Supporting Code

for “Animal-vectored nutrient flows over resource gradients influence the nature of local and meta-ecosystem functioning”

Introduction

This notebook contains the code used to run numerical analyses of a meta-ecosystem model in which consumers move across ecosystem borders in different ways with respect to resource gradients.

In section Model Numerical Analyses, we perform numerical analyses of our model of consumer movement in a two-patch meta-ecosystem, using the set of equilibrium equations derived in Mathematica (see Mathematica notebook in this same repository; Wolfram Research Inc. (2020)). In section Changes to biomass distribution in local and meta-ecosystem, we investigate how biomass distribution changes following different types of consumer movement and how this may affect top-down influences of consumer movement on local and meta-ecosystem functions. Sections Visual Deliverables and Summary Tables contain code to reproduce the figures and tables shown in the main text and Supplementary Materials.

In this introduction, we briefly define the contrasting consumer movement types we investigate with this model (Defining movement types), discuss separation of scales in time and space (Separation of Spatial and Temporal scales), describe the state variables and parameters in the model (State variables and parameters of the model), and describe the analyses setup (Analyses setup).

Defining movement types

Here, we focus on consumers moving over gradients of inorganic nutrients availability between two local ecosystems. Consumer movement—either gradient neutral, along-, or against-gradient—connects these two local ecosystems in a meta-ecosystem (Loreau, Mouquet, and Holt 2003). We define consumer movement as gradient neutral when it happens among ecosystems that have equal nutrient availability. Along-gradient movement takes place when consumers move following the gradient direction. Conversely, against-gradient consumer movement runs counter to the direction of the gradient. Consistent with recent reviews (Gounand et al. 2018), in the following sections we also use diffusive and non-diffusive to refer to along- and against-gradient consumer movement, respectively. In our model, we create resource gradients by varying the relative values of the inorganic nutrient input rates in the two ecosystems comprising our meta-ecosystem (parameter I_i; see Table 1 in the main text or below for a list of model parameters and state variables).

Separation of spatial and temporal scales

While organisms can at times cross into an ecosystem that is equally suitable for them immediately after leaving their original one, this is not often the case. Rather, organisms spend significant amounts of time searching for a new suitable ecosystem, traversing the so-called “unsuitable” matrix (Weisser and Hassell 1996). Models of organismal movement in meta-ecosystems should account for this time spent outside ecosystem borders, as it may lead to loss of biomass and nutrients from the system (Gounand et al. 2018).

We model consumer movement in our meta-ecosystem model by combining a dispersers’ pool (Weisser and Hassell 1996) with the time scales separation (TSS) mathematical technique (Otto and Day 2011). Briefly, adding a dispersers’ pool to a classic meta-ecosystem model means that consumer do not instantly travel from donor to recipient ecosystem. Instead, they move from the donor to the dispersers’ pool and then from there to the recipient (see Figure 1c in the main text for a visual representation). This allows us to separate ecosystem processes that happen locally from those that happen regionally.

Some of these ecosystem processes, like primary productivity in terrestrial ecosystems, take place over long time frames (slow processes). Others, like consumer movement, happen over short time frames (fast processes). TSS explicitly accounts for these differences by allowing fast processes to reach a quasi-equilibrium while holding slow ones constant. The quasi-equilibrium is then plugged into the equation(s) for the slow processes to investigate how these processes evolve as the fast ones stay constant (Otto and Day 2011).

The combined use of the dispersers’ pool and TSS allows us to separate local and regional dynamics in our system and to investigate the effects of consumer movement at either spatial extent at the same time. See section Model Derivation in the main text for further details on our novel approach to model consumer movement in meta-ecosystems.

State variables and parameters of the model

The table below lists the model’s state variables and parameters, and is analogous to Table 1 in the main text. For a visual representation of the model, please see Figure 1c in the main text.

Variable	Description	Units	Range
Ni	Nutrients stocks in patch i	g	>0
Pi	Producers stocks in patch i	g	>0
Ci	Consumers stocks in patch i	g	>0
Q	Dispersers’ Pool	g	>0

Parameter	Description	Units	Range
Ii	Inorganic nutrient input rate into patch i	g*t-1	[0,20]
\(l\)	Inorganic nutrient output rate	t-1	[0,10]
ui	Producer uptake rate in patch i	(g*t)-1	[0,10]
ai	Consumer attack rate in patch i	(g*t)-1	[0,10]
\(\epsilon\)i	Consumer assimilation efficiency in patch i	unitless	[0,1]
hi	Biomass loss rate from Pi	t-1	[0,10]
di	Biomass loss rate from Ci	t-1	[0,10]
g	Movement rate from C1 to Q	t-1	[0,10]
m	Movement rate from Q to C2	t-1	[0,10]
c	Biomass loss rate from Q	t-1	[0,10]

We derived the model’s single feasible equilibrium in Mathematica (Wolfram Research Inc. (2020); see Model Derivation in the main text for further details).

Analyses setup

This notebook assumes the following folder structure:

.
├── Code
├── Data
└── Results

This structure allows us to take advantage of relative paths when loading data or saving figures. Note that . stands for the root folder. If you wish to recompile this notebook, please organize your folder structure accordingly.

This notebook uses the following packages:

Show code

library(renv) # project dependency (packages, mostly) management
library(tidyverse) # collection of packages improving base R 
library(plotrix) # additional plotting functions
library(lhs) # use latin hypercube to create toy datasets
library(ggpubr) # extra themes for ggplot2
library(ggpol) # extra geoms for ggplot2
library(patchwork) # easily create and layout multi-panel plots
library(gt) # building tables in a tidyverse framework
library(scales) # scaling functions for ggplot2
library(khroma) # color-blind palettes for ggplot2
library(xaringanExtra) # extra tabbing functionality for rmarkdown documents
library(tictoc) # timing 
library(RColorBrewer) # color palettes manipulation
library(MetBrewer) # additional color palettes
library(broom) # data wrangling

# Packages that need to be installed but do not require loading

# library(distill) # notebook templates
# library(rmarkdown) # build interactive documents in R 

In the interest of reproducibility, we use package renv to keep track of this notebook’s packages.

Show code

renv::snapshot()

* Lockfile written to '~/PhD/Chapter_3/Code/renv.lock'.

We compiled this notebook on a machine using the following version of R, operative system, and necessary packages:

Show code

utils:::print.sessionInfo(sessionInfo()[-8])

R version 4.2.0 (2022-04-22)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.4

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
 [1] broom_0.8.0         MetBrewer_0.2.0     RColorBrewer_1.1-3 
 [4] tictoc_1.0.1        xaringanExtra_0.6.0 khroma_1.8.0       
 [7] scales_1.2.0        gt_0.6.0            patchwork_1.1.1    
[10] ggpol_0.0.7         ggpubr_0.4.0        lhs_1.1.5          
[13] plotrix_3.8-2       forcats_0.5.1       stringr_1.4.0      
[16] dplyr_1.0.9         purrr_0.3.4         readr_2.1.2        
[19] tidyr_1.2.0         tibble_3.1.7        ggplot2_3.3.6      
[22] tidyverse_1.3.1     renv_0.15.5        

Model numerical analyses

Here, we perform the analyses of the meta-ecosystem model described in section Numerical Analyses of the main text. If you are interested in the code used to produce the main text’s figures and tables, this can be found in sections Visual Deliverables and Summary Tables below.

We conduct three rounds of analyses, one each for the three types of consumer movement we are interested in: gradient-neutral, along-, and against-gradient. Each round of analysis comprises iterating the model over n = 10 000 parameter sets, running stability analyses, and graphing the results. The analyses will follow these steps:

1. Fit randomly-drawn parameter sets to the model’s equilibrium and solve for state variable values
2. Calculate meta-ecosystem function values
3. Flag parameter sets that produce equilibria < 0, hence without biological meaning
4. Run stability analyses on each parameter set to identify those producing unstable results
5. Exclude parameter sets that are unstable and/or have no biological meaning from further analyses
6. Explore results graphically

In the code chunk below, we generate two sets of randomly-drawn parameter values to use in all scenarios. We use function lhs() to do so. One set includes parameters values scaled 0–10 (DATA1), while the other one includes parameter values scaled 0–1 (DATA2).

Show code

# Below we generate random values for those parameters that will vary
# during the simulations: attack rates, efficiency, death rates, and movement 
# rates

DATA1 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for(i in seq(1,100,1)){
  
  # Use lhs function in R. 100 bins, 11 iterations - one for each parameter that
  # is scaled between 0 and 10.
  loop_1 <- randomLHS(100, 11)
  
  # Use rbind to stack each iteration on top of the other to create a one column
  # list of 10000 samples with equal distribution in 100 bins.
  # Multiplied by 10 to get range of data from 0 to 10
  DATA1 <- rbind(DATA1, data.frame(i=i, Result=loop_1*10))
  
  # print the i value so you can track progress
  # print(i)
  # Close the loop
}

# This is for 2 parameters that is scaled 0 to 1.
DATA2 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for(i in seq(1,100,1)){

  # Use lhs function in R. 100 bins, 2 iterations.
  loop_2 <- randomLHS(100, 2)
  
  # Use rbind to stack each iteration on top of the other to create a one column 
  # list of 200000 samples with equal distribution in 100 bins.
  # Multiplied by 10 to get range of data from 0 to 10
  DATA2 <- rbind(DATA2, data.frame(i=i, Result=loop_2))
  
  # print the i value so you can track progress
  # print(i)
  # Close the loop
}

# Save these two datasets for later re-use
# write.csv(DATA1, file = "../Data/DATA1.csv", row.names = FALSE)
# write.csv(DATA2, file = "../Data/DATA2.csv", row.names = FALSE)

Alternatively, we can load the two datasets provided alongside this R notebook, in folder ../Data/. These were similarly generated using function lhs(), scaled either 0–10 or 0–1, and then used to produce the results and figures presented in the paper.

Loading the provided datasets is the default approach. To generate new randomly drawn parameter values, set eval=FALSE in the chunk below and eval=TRUE in the code chunk above.

Show code

DATA1 <- read.csv(file = "../Data/DATA1.csv", header = TRUE)

DATA2 <- read.csv(file = "../Data/DATA2.csv", header = TRUE)

Control scenario: gradient-neutral consumer movement

We begin our numerical analyses of our model with the control scenario, in which consumers move between two ecosystems that have equal environmental fertility. Consumer movement in this scenario is gradient-neutral, as the environment is homogeneous in terms of resource availability. The results of this first set of model iterations produce the values used at the denominator of the response ratio (LRR) used to capture changes in local and meta-ecosystem function values (see section Visual Deliverables below). Below, we describe each step of the analyses in detail: later analyses will not describe these steps again.

Step 1: plug the parameter sets into the model

First, we allocate an empty dataframe PRED to contain the results of each iteration of the model. This dataframe will contain the parameter values and the resulting local and meta-ecosystem function values.

Show code

# Allocate an empty data frame to save the data in, then allocate some rows and 
# columns. We will replace the current numbers with results below 
# Data we are calculating are (bio)mass for soil nutrients (N1, N2), producers 
# (P1, P2), consumers (C1, C2), disperses (Q)

PRED <- NULL

PRED <-rbind(PRED, data.frame(TIME=seq(0,100,1), I1 = seq(0,100,1), I2 = seq(0,100,1), 
                              l = seq(0,100,1), u1 = seq(0,100,1), u2 = seq(0,100,1), 
                              a1 = seq(0,100,1), a2 = seq(0,100,1), h1 = seq(0,100,1), 
                              h2 = seq(0,100,1), d1 = seq(0,100,1), d2 = seq(0,100,1), 
                              g = seq(0,100,1), m = seq(0,100,1), c = seq(0,100,1), 
                              e1 = seq(0,100,1), e2 = seq(0,100,1), N1 = seq(0,100,1), 
                              P1 = seq(0,100,1), C1 = seq(0,100,1), N2 = seq(0,100,1), 
                              P2 = seq(0,100,1), C2 = seq(0,100,1), Q = seq(0,100,1), 
                              FLUX_P1 = seq(0,100,1), FLUX_C1 = seq(0,100,1), 
                              FLUX_P2 = seq(0,100,1), FLUX_C2 = seq(0,100,1), 
                              FLUX_Eco_1 = seq(0,100,1), FLUX_Eco_2 = seq(0,100,1), 
                              FLUX_Q = seq(0,100,1), PROD_P1 = seq(0,100,1), 
                              PROD_C1 = seq(0,100,1), PROD_P2 = seq(0,100,1), 
                              PROD_C2 = seq(0,100,1)))

Now we run the model. We first set the values for both nutrient input rates in ecosystem 1 (I1 in the code chunks) and ecosystem 2 (I2) at 10. Recall that this is the control scenario with equal environmental fertility. We also set the value of the inorganic nutrient leaching rate (l) at 0.1.

Note that we refer to I1 and I2 as either the nutrient input rates or the environmental fertility of the two ecosystems

The code chunk below runs the model using a for loop that, over 10 000 time steps, iterates the model’s single feasible equilibrium to calculate:

• stocks of Nutrients (N1, N2), Producers (P1, P2), and Consumers (C1, C2)
• the Dispersers’ Pool (Q) stock
• the recycling flux for each biotic compartment (FLUX_X, where X is the compartment)
• the productivity for each biotic compartment (PROD_X, where X is the trophic compartment)

in both local ecosystems. At each time step, we draw a random value for each parameter from either the DATA1 or DATA2 dataset—depending on whether the parameter is scaled 0–10 or 0–1, respectively. Then, we plug these values into the model’s equilibria and ecosystem function formulas derived in Mathematica (Wolfram Research Inc. 2020), and finally store both the randomly drawn parameter values and ecosystem function values in PRED.

The feasibility conditions for this model are too complex to solve analytically. Hence, we are not checking that the randomly-drawn parameters values satisfy them. See below for stability analyses of each parameter set.

Show code

# Here we set values for parameters that we will not vary during the simulations
# Nutrients input rate in each patch
I1 = 10
I2 = 10
# Patch leaching rate/soil loss rate
l = 0.1

# We want to start simulations at 0
i1 = 0

for (i in seq(0,9999,1)){
  # go to the next row
  i1 = i1 + 1

# Let's sample the parameter values from DATA. We will sample these i times as 
# defined by the for loop above

  # producers uptake rates
  u1 = DATA1[i1,2] # in ecosystem 1
  u2 = DATA1[i1,3] # in ecosystem 2
  # consumers attack rates
  a1 = DATA1[i1,4] # in ecosystem 1
  a2 = DATA1[i1,5] # in ecosystem 2
  # death rates (= recycling rate)
  h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
  h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystems 2
  d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
  d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
  c = DATA1[i1,10] # death rate in the disperser's pool Q
  # movement rate to the Dispersers' pool Q
  g = DATA1[i1,11]
  # movement rate from the Dispersers' pool Q
  m = DATA1[i1,12]
  # assimilation efficiency of consumers 
  e1 = DATA2[i1,2] # in ecosystem 1
  e2 = DATA2[i1,3] # in ecosystem 2
  
  PRED[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
                 # Nutrient stock in ecosystem 1
                 ((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Primary Producers biomass in ecosystem 1
                 (d1 + g)/(a1*e1),
                 # Consumers biomass in ecosystem 1
                 (e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Nutrient stock in ecosystem 2
                 (-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
                 # Primary Producers biomass in ecosystem 2
                 (l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 # Consumer biomass in ecosystem 2
                 ((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
                 # Consumer biomass in the Dispersers' Pool (Q)
                 (e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 # Nutrient Flux for Primary Producers in Ecosystem 1
                 ((d1 + g)*h1)/(a1*e1),
                 # Nutrient Flux for Consumers in Ecosystem 1
                 (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Nutrient Flux for Primary Producers in Ecosystem 2
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 # Nutrient Flux for Consumers in Ecosystem 2
                 (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 # Nutrient Flux for Ecosystem 1
                 ((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Nutrient Flux for Ecosystem 2
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 # Nutrient Flux in the Dispersers' Pool (Q)
                 (c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 # Primary Productivity in Ecosystem 1
                 ((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 # Consumer (=Secondary) Productivity in Ecosystem 1
                 (e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 # Primary Productivity in Ecosystem 2
                 ((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)),
                 # Consumer (=Secondary) Productivity in Ecosystem 2
                 (e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
                 )
}

Step 2: calculate meta-ecosystem functions values

To investigate any influence of consumer movement on nutrient cycling at regional spatial scales, we calculate local and meta-ecosystem total biomass (Ntot,Ptot, Ctot), meta-ecosystem nutrient flux (MetaEcoFlux) and trophic compartment nutrient flux (FLUX_Ptot, FLUX_Ctot) and productivity (PROD_Ptot, PROD_Ctot). Furthermore, we manually calculate each ecosystem’s own total recycling flux as the sum of the flux of its 2 biotic compartments, as a check on the Mathematica-derived formulas used above.

Show code

# calculate meta-ecosystem biomass stocks for each trophic compartment,
# a.k.a., meta-ecosystem biomass
PRED$Ntot <- PRED$N1+PRED$N2
PRED$Ptot <- PRED$P1+PRED$P2
PRED$Ctot <- PRED$C1+PRED$C2

# now calculate recycling flux for the local ecosystem manually to double check 
# the formulae used above are correct 
PRED$FLUX_Eco1_check <- PRED$FLUX_P1 + PRED$FLUX_C1
PRED$FLUX_Eco2_check <- PRED$FLUX_P2 + PRED$FLUX_C2

# and the meta-ecosystem recycling flux
PRED$FLUX_Ptot <- PRED$FLUX_P1 + PRED$FLUX_P2
PRED$FLUX_Ctot <- PRED$FLUX_C1 + PRED$FLUX_C2
PRED$MetaEcoFlux <- PRED$FLUX_Eco_1 + PRED$FLUX_Eco_2 - PRED$FLUX_Q

# finally, calculate the production for each compartment at meta-ecosystem level
PRED <- dplyr::mutate(PRED, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PRED <- dplyr::mutate(PRED, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")

Step 3: flag parameter sets with no biological meaning

Parameter sets that produce values of the model’s state variable that are \(\leq 0\) lack biological meaning. Accordingly, in the code chunk below, we flag these parameter sets for later exclusion from the analysis.

Show code

# Equilibrium stocks less than or equal to 0 make no biological sense, so we 
# will flag them in the PRED dataset.
PRED <- PRED %>% dplyr::mutate(., biosense = ifelse(N1 > 0 & 
                                                    P1 > 0 & 
                                                    C1 > 0 & 
                                                    N2 > 0 & 
                                                    P2 > 0 & 
                                                    C2 > 0 & 
                                                    Q > 0, 
                                                    "yes", 
                                                    "no"), 
                               .after = MetaEcoFlux) %>% 
  dplyr::mutate_at(vars(biosense), factor)

Step 4: stability analyses

We perform the stability analyses on the model by fitting the same sets of parameter values and equilibrium state variable values used above to the partial derivatives extracted from the model’s Jacobian matrix evaluated in Mathematica (Wolfram Research Inc. 2020). We evaluate stability of each parameter set by checking whether the real part of the leading eigenvalue of the Jacobian is positive or negative (Otto and Day 2011). If the real part of the Jacobian’s leading eigenvalue is negative, the equilibrium is stable (Otto and Day 2011).

Show code

# create an empty dataframe
MathStab <- NULL

for (i in 1:nrow(PRED)) {
  
  # fetch parameter values from the equilibrium simulations df
  I1 = PRED$I1[i]
  I2 = PRED$I2[i]
  l = PRED$l[i]
  u1 = PRED$u1[i]
  u2 = PRED$u2[i]
  a1 = PRED$a1[i]
  a2 = PRED$a2[i]
  h1 = PRED$h1[i]
  h2 = PRED$h2[i]
  d1 = PRED$d1[i]
  d2 = PRED$d2[i]
  g = PRED$g[i]
  m = PRED$m[i]
  c = PRED$c[i]
  e1 = PRED$e1[i]
  e2 = PRED$e2[i]
  
  # Put in the solutions to the equilibrium values here:
  dN1 = PRED$N1[i]
  dP1 = PRED$P1[i]
  dC1 = PRED$C1[i]
  dN2 = PRED$N2[i]
  dP2 = PRED$P2[i]
  dC2 = PRED$C2[i]
  
  # create the Jacobian from Mathematica
  # NOTE: in transposing from Mathematica, the Jacobian is here assembled by row
  Jacob <- rbind(c(-l-dP1*u1, h1-dN1*u1, d1, 0, 0, 0), 
                 c(dP1*u1, -a1*dC1-h1+dN1*u1, -a1*dP1, 0, 0, 0), 
                 c(0, a1*dC1*e1, -d1-g+a1*e1*dP1, 0, 0, 0), 
                 c(0, 0, 0, -l-dP2*u2, h2-dN2*u2, d2), 
                 c(0, 0, 0, dP2*u2, -a2*dC2-h2+dN2*u2, -a2*dP2), 
                 c(0, 0, (g*m)/(c+m), 0, a2*dC2*e2, -d2+a2*e2*dP2))
  
  # Check to see if the real part of the leading eigenvalue is less than 0 or not - if it is than the system is stable.
  if (max(Re(base::eigen(Jacob)$values)) < 0 ){ 
    stable <- "stable"
  } else {
    stable <- "unstable"
  }
  
  MathStab <- rbind(MathStab, data.frame(TIME=i, I1, I2, l, u1, u2, a1, a2, h1,
                                         h2, d1, d2, g, m, c, e1, e2, dN1, dP1,
                                         dC1, dN2, dP2, dC2,
                                         EiV1 = Re(eigen(Jacob)$values[1]),
                                         EiV2 = Re(eigen(Jacob)$values[2]),
                                         EiV3 = Re(eigen(Jacob)$values[3]),
                                         EiV4 = Re(eigen(Jacob)$values[4]),
                                         EiV5 = Re(eigen(Jacob)$values[5]),
                                         EiV6 = Re(eigen(Jacob)$values[6]), 
                                         maxEv = max(Re(base::eigen(Jacob)$values)),
                                         stable = stable,
                                         biosense = PRED$biosense[i]))
}

MathStab$stable <- as.factor(MathStab$stable)

# separate unstable equilibria to work with later
MathStabUS <- subset(MathStab, MathStab$stable == "unstable")

According to the stability analyses run above, about 1.55% of the parameter sets produced unstable results (i.e., 155 out of 10000 iterations).

Step 5: exclude unstable parameter sets from further analyses

Now, let’s check if the parameter sets that produce unstable equilibria are also those that produce equilibria lacking biological meaning—i.e., those where state variables values at equilibrium are \(\leq 0\).

Show code

biononsense <- subset(PRED, PRED$biosense == "no")

identical(as.numeric(MathStabUS[ , "TIME"]), biononsense[ , "TIME"])

[1] TRUE

It appears that they are. Hence, let’s add the stable/unstable information to PRED.

Show code

PRED <- left_join(PRED, select(MathStab, !c(I1:maxEv)), by = c("TIME", "biosense"))

Now, let’s exclude these unstable parameter sets from the analyses below and from our results.

We are also going to sample the dataset PRED, that contains the results of our model’s iterations, for 1000 random iterations results. We will store these in PRED_1k and use them in Figures 1, 2, 3, 4, 5 below to produce some preliminary graphs of our results.

Show code

PREDpos <- filter(PRED, PRED$biosense == "yes" & PRED$stable == "stable")

# sample PREDpos to only use 1000 random simulation results
PRED_sample1000 <- PREDpos[sample(nrow(PREDpos), size = 1000), ]

PRED_1k <- droplevels(PRED_sample1000)

Before proceeding further, let’s save the PRED object to disk, for ease of use in the future.

Show code

write.csv(PRED, file = "../Results/PRED.csv", row.names = FALSE)

Step 6: explore results graphically

We now investigate how the presence of unidirectional, gradient-neutral consumer movement influences the meta-ecosystem’s behavior using box-plots produced using PRED_1k. Note that in the figures below, “Donor” always indicates ecosystem 1 and “Recipient” always indicates ecosystem 2.

Biomass and nutrient stock accumulation

Show code

# pivot the dataset to longer format for boxplot plotting
PRED_biomass_long <- select(PRED_1k, N1:Q, Ntot:Ctot) %>%
  pivot_longer(.,
               names_to = "Compartment",
               values_to = "Stock",
               cols = c(N1, P1, C1, N2, P2, C2, Q,
                        Ntot, Ptot, Ctot)) %>%
  dplyr::mutate(.,
                Ecosystem = if_else(Compartment == "N1" |
                                      Compartment == "C1" |
                                      Compartment == "P1",
                                    "Donor",
                                    if_else(Compartment == "N2" |
                                              Compartment == "P2" |
                                              Compartment=="C2",
                                            "Recipient",
                                            if_else(Compartment == "Q",
                                                    "Dispersers Pool",
                                                    "Meta-ecosystem"))))

PRED_biomass_long$Compartment <- as_factor(PRED_biomass_long$Compartment)
PRED_biomass_long$Ecosystem <- as_factor(PRED_biomass_long$Ecosystem)

comparts_medians <- PRED_biomass_long %>%
  dplyr::group_by(Compartment, Ecosystem) %>%
  dplyr::summarise(., median = median(Stock),
                   .groups = "keep")

ggplot(data = PRED_biomass_long,
       aes(x = Compartment, y = Stock, fill = Ecosystem, col = Ecosystem)) +
  geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE,
                 outlier.shape = 25, jitter.shape = 21) +
  geom_label(data = comparts_medians,
             aes(y = median, label = round(median,2)),
             col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
             nudge_x = -0.18, label.padding = unit(0.15, "lines"),
             label.size = 0.15) +
  facet_grid(.~Ecosystem, scales = "free") +
  scale_color_manual(values = met.brewer("Egypt", 4)) +
  scale_fill_manual(values = met.brewer("Egypt", 4)) +
  theme_pubr() +
  coord_cartesian(y = c(0, 100)) +
  theme(legend.position = "none")

Figure 1: Biomass stock values at equilibrium for local and meta-ecosystem scale, drawn using 1k randomly-drawn data points from the model simulations. Consumer movement from ecosystem 1 to ecosystem 2 via the dispersers’ pool results in foraging pressure of consumers being transferred from producers in ecosystem 1 to those in ecosystem 2. Primary producers in ecosystem 1 are thus released from foraging pressure while those in ecosystem 2 are depressed. In turn, this leads to differences in nutrient stocks between the two ecosystems at both the local and meta-ecosystem extents. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations (i.e., 10k iterations). Values on the y-axis are limited between [0, 100] to zoom in on the variation in the results.

The graph above shows a spatial trophic cascade (Knight et al. 2005). In ecosystem 2, the influx of consumers from ecosystem 1 increases the consumer biomass, in turn depressing the producers’ abundance and releasing the nutrient stock from the producers’ trophic pressure. Conversely, in ecosystem 1, consumers’ numbers fall steadily—due to the combined action of mortality and emigration towards ecosystem 2—thus releasing the producers from the trophic pressure exerted by consumers. In turn, this lets producers grow unchecked, depressing the nutrient stocks in ecosystem 1. We can see this as well by looking at the change of log₁₀ biomass values for each trophic compartment with respect to the two movement parameters: g, the rate of consumer movement from ecosystem 1 to the dispersers’ pool Q (Figure 2), and m, the rate of consumer movement from the dispersers’ pool Q to ecosystem 2 (Figure 3).

Show code

eco_levels <- c("N1" = "Nutrient", "P1" = "Producers", "C1" = "Consumers",
                "N2" = "Nutrient", "P2" = "Producers", "C2" = "Consumers",
                "Q" = "Dispersers", "Ntot" = "Nutrient",
                "Ptot" = "Producers", "Ctot" = "Consumers")

# pivot the dataset to longer format for boxplot plotting
PRED_long <- pivot_longer(PRED_1k,
                          names_to = "Compartment",
                          values_to = "Stock",
                          cols = c(N1, P1, C1, N2, P2, C2, Q,
                                   Ntot, Ptot, Ctot)) %>%
  dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" |
                                         Compartment == "C1" |
                                         Compartment == "P1",
                                       "Donor",
                                       if_else(Compartment == "N2" |
                                                 Compartment == "P2" |
                                                 Compartment=="C2", "Recipient",
                                               if_else(Compartment == "Q",
                                                       "Dispersers",
                                                       "Meta-ecosystem"))))

PRED_long$Compartment <- as_factor(PRED_long$Compartment)
PRED_long$Ecosystem <- as_factor(PRED_long$Ecosystem)

ggplot(PRED_long, aes(g, log10(Stock), col = Ecosystem)) +
  geom_point(alpha = 0.15) +
  geom_smooth(method = "lm", col = "white", alpha = 1, se = T, lwd = .4) +
  facet_grid(.~Compartment) +
  scale_color_manual(values = met.brewer("Egypt", 4)) +
  theme_pubr() +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
        text = element_text(size = 11))

Figure 2: Change in biomass in the meta-ecosystems’ trophic compartments with change in the movement rate from the donor ecosystem 1 to the dispersers’ pool (g). Each facet presents results for a trophic compartment, with color identifying local and meta-ecosystem. White lines are regression lines with 95% CI (grey shaded areas). Note the log₁₀ scale on the y-axis. Graph produced using 1000 randomly-drawn data points from the model’s analyses.

Show code

ggplot(PRED_long, aes(m, log10(Stock), col = Ecosystem)) +
  geom_point(alpha = 0.15) +
  geom_smooth(method = "lm", col = "white", alpha = 1, se = T, lwd = .4) +
  facet_grid(.~Compartment) +
  scale_color_manual(values = met.brewer("Egypt", 4)) +
  theme_pubr() +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
        text = element_text(size = 11))

Figure 3: Change in biomass in the meta-ecostystem’s trophic compartments with change in movement rate from the disperser pool to the recipient ecosystem 2 (m). Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log₁₀ scale on the y-axis. All other specifications as in Figure 2.

Finally, we can see a similar effect when looking at the change in biomass in the two local ecosystems with change in the mortality rate of consumers while in the dispersers’ pool (parameter c in the model; Figure 4).

Show code

ggplot(PRED_long, aes(c, log10(Stock), col = Ecosystem)) +
  geom_point(alpha = 0.15) +
  geom_smooth(method = "lm", col = "white", alpha = 1, se = T, lwd = .4) +
  facet_grid(.~Compartment) +
  scale_color_manual(values = met.brewer("Egypt", 4)) +
  theme_pubr() +
  guides(colour = guide_legend(override.aes = list(alpha = 1))) +
  theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
        text = element_text(size = 11))

Figure 4: Change in biomass in the meta-ecostystem’s trophic compartments with change in the consumer death rate while in the dispersers’ pool (c). This could represent predation risk, but also environmental hazards and stochastic events such as crossing a dangerous linear feature of the landscape—e.g., a busy highway. Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log₁₀ scale on the y-axis. All other specifications as in Figure 2.

Nutrient flux

In the graphs below (Figure 5), we look at nutrient flux values in each ecosystem and in the meta-ecosystem for each iteration of the model.

First, we will extract the flux columns (see above) from the PRED_1k dataframe into a new one which we will call FluxPRED. We will then pivot this dataframe from wide to long format, using function pivot_longer from package tidyr.

Show code

FluxPRED <- PRED_1k %>%
  dplyr::select(TIME:c,FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
  tidyr::pivot_longer(., names_to = "Compartment",
                      values_to = "Flux",
                      cols = c(FLUX_P1:MetaEcoFlux)) %>%
  dplyr::mutate(., Scale = if_else(Compartment == "FLUX_P1" |
                                     Compartment == "FLUX_C1" |
                                     Compartment == "FLUX_Eco_1",
                                   "Donor",
                                   if_else(Compartment == "FLUX_P2" |
                                             Compartment == "FLUX_C2" |
                                             Compartment == "FLUX_Eco_2",
                                           "Recipient",
                                           "Meta-ecosystem"))) %>%
  dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))

Now, we produce the box-and-whiskers plot of recycling flux in local and meta-ecosystem.

Show code

FluxPRED$Scale <- as.factor(FluxPRED$Scale)
FluxPRED$Compartment <- as.factor(FluxPRED$Compartment)

comparts_medians <- FluxPRED %>%
  dplyr::group_by(Compartment, Scale) %>%
  dplyr::summarise(., median = median(Flux), .groups = "keep")

FluxPRED %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment,
                                             c("FLUX_P1", "FLUX_C1", "FLUX_Eco_1",
                                               "FLUX_P2", "FLUX_C2", "FLUX_Eco_2",
                                               "MetaEcoFlux"))) %>%
  dplyr::mutate(., Scale = fct_relevel(Scale,
                                       c("Donor",
                                         "Recipient",
                                         "Meta-ecosystem"))) %>%
  ggplot(., aes(x = Compartment, y = Flux)) +
  geom_boxjitter(outlier.intersect = TRUE, alpha = 0.25, outlier.shape = 25, jitter.shape = 21) +
  geom_label(data = comparts_medians, aes(y = median, label = round(median,2)),
             col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
             nudge_x = -0.18, label.padding = unit(0.15, "lines"),
             label.size = 0.15) +
  facet_grid(.~Scale, scales = "free") +
  scale_x_discrete(labels = c("FLUX_C1" = "C1", "FLUX_P1" = "P1",
                              "FLUX_Ptot" = "Producers", "FLUX_C2" = "C2",
                              "FLUX_P2" = "P2", "FLUX_Ctot" = "Consumers",
                              "FLUX_Eco_1" = "Donor", "FLUX_Eco_2" = "Recipient",
                              "MetaEcoFlux" = "Meta-ecosystem")) +
  theme_pubr() +
  coord_cartesian(y = c(0, 500)) +
  theme(legend.position = "none", axis.title.x = element_blank())

Figure 5: Nutrient flux at local and meta-ecosystem scales when consumers move from ecosystem 1 to 2. The two ecosystems have the same fertility—i.e., basal inorganic input values are the same. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations. The y-axis is constrained between [0, 500] to show the variation in the results. Graph produced using 1000 randomly-drawn data points from the model’s analyses.

Experimental scenarios: along- and against-gradient consumer movement

Here, we investigate the influence of consumers moving from ecosystem 1 to ecosystem 2 when there is a gradient of nutrient availability between the two ecosystems. We will investigate two scenarios: higher environmental fertility in the donor ecosystem 1 (i.e., I1 >> I2) and higher environmental fertility in the recipient ecosystem 2 (i.e., I1 << I2). Note that the environmental leaching rate of each ecosystem—i.e., the rate at which nutrients leave local ecosystems independently of trophic processes—is invariant and equal between the two ecosystems (l = 0.1).

The analysis will follow the same steps as before: draw parameter values from our DATA1 and DATA2 data objects, simulate model behavior, run stability analyses, remove equilibria that are unstable and/or lack biological sense, and graphically inspect and interpret the results.

Along-gradient consumer movement

We begin by simulating a meta-ecosystem where the donor ecosystem 1 has higher environmental fertility than the recipient ecosystem 2. Consumers move along-gradient, in a diffusive way that resembles how organismal movement has been represented in previous studies (e.g., Gravel et al. 2010).

Show code

# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 > I2
I1=18
I2=2

# Leaching does not change as is still equal across ecosystems
l = 0.1

# Allocate an empty data frame to save the data in
PREDI2 <- NULL
# Allocate some rows and columns. We will replace the current numbers with results below
# Data we are calculating are (bio)mass for soil nutrients (N1, N2), producers (P1, P2), consumers (C1,C2), disperses (Q)

PREDI2 <-rbind(PREDI2, data.frame(TIME=seq(0,100,1), 
                                  I1 = I1,  I2 = I2, l = l, 
                                  u1 = seq(0,100,1), u2 = seq(0,100,1), 
                                  a1 = seq(0,100,1), a2 = seq(0,100,1), 
                                  h1 = seq(0,100,1), h2 = seq(0,100,1), 
                                  d1 = seq(0,100,1), d2 = seq(0,100,1), 
                                  g = seq(0,100,1), m = seq(0,100,1), 
                                  c = seq(0,100,1), 
                                  e1 = seq(0,100,1), e2 = seq(0,100,1), 
                                  N1=seq(0,100,1), P1=seq(0,100,1), 
                                  C1=seq(0,100,1), N2=seq(0,100,1), 
                                  P2=seq(0,100,1), C2=seq(0,100,1), 
                                  Q = seq(0,100,1), 
                                  FLUX_P1 = seq(0,100,1), 
                                  FLUX_C1 = seq(0,100,1), 
                                  FLUX_P2 = seq(0,100,1), 
                                  FLUX_C2 = seq(0,100,1), 
                                  FLUX_Eco_1 = seq(0,100,1), 
                                  FLUX_Eco_2 = seq(0,100,1), 
                                  FLUX_Q = seq(0,100,1), 
                                  PROD_P1 = seq(0,100,1), 
                                  PROD_C1 = seq(0,100,1), 
                                  PROD_P2 = seq(0,100,1), 
                                  PROD_C2 = seq(0,100,1)))

# We want to start simulations at 0
i1 = 0

for (i in seq(0,9999,1)){
  # go to the next row
  i1 = i1 + 1

  # Let's sample the parameter values from DATA. We will sample these x times as defined by the i for loop above
  # producers uptake rates
  u1 = DATA1[i1,2] # in ecosystem 1
  u2 = DATA1[i1,3] # in ecosystem 2
  # consumers attack rates
  a1 = DATA1[i1,4]
  a2 = DATA1[i1,5]
  # death rates (= recycling rate)
  h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
  h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystem 2
  d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
  d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
  c = DATA1[i1,10] # death rate in the disperser's pool Q
  # movement rate to the Disperser's pool Q
  g = DATA1[i1,11]
  # movement rate from the Disperser's pool Q
  m = DATA1[i1,12]
  # efficiency of consumers 
  e1 = DATA2[i1,2] # in ecosystem 1
  e2 = DATA2[i1,3] # in ecosystem 2
  
  PREDI2[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
                 ((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (d1 + g)/(a1*e1),
                 (e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
                 (l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 ((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
                 (e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 ((d1 + g)*h1)/(a1*e1),
                 (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 ((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 (c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 ((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 (e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 ((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l -I1*u1)))*u2)),
                 (e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
                 )
  }

# calculate meta-ecosystem biomass stocks for each trophic compartment,
# a.k.a., meta-ecosystem productivity
PREDI2$Ntot <- PREDI2$N1 + PREDI2$N2
PREDI2$Ptot <- PREDI2$P1 + PREDI2$P2
PREDI2$Ctot <- PREDI2$C1 + PREDI2$C2

# now calculate recycling flux for the local ecosystem
PREDI2$FLUX_Eco1_check <- PREDI2$FLUX_P1 + PREDI2$FLUX_C1
PREDI2$FLUX_Eco2_check <- PREDI2$FLUX_P2 + PREDI2$FLUX_C2

# and the meta-ecosystem recycling flux
PREDI2$FLUX_Ptot <- PREDI2$FLUX_P1 + PREDI2$FLUX_P2
PREDI2$FLUX_Ctot <- PREDI2$FLUX_C1 + PREDI2$FLUX_C2
PREDI2$MetaEcoFlux <- PREDI2$FLUX_Eco_1 + PREDI2$FLUX_Eco_2 - PREDI2$FLUX_Q

# Finally, calculate the production for each compartment at the meta-ecosystem scale
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")

# Equilibrium stocks less than or equal to 0 make no biological sense, so we 
# will flag them in the PREDI2 dataset.
PREDI2 <- PREDI2 %>% 
  dplyr::mutate(., biosense = ifelse(N1 > 0 & 
                                     P1 > 0 & 
                                     C1 > 0 & 
                                     N2 > 0 & 
                                     P2 > 0 & 
                                     C2 > 0 & 
                                     Q > 0, 
                                     "yes", 
                                     "no"), 
                .after = MetaEcoFlux) %>% 
  dplyr::mutate_at(vars(biosense), factor)

Stability Analyses

Here we perform stability analyses for the model in this scenario where I1 >> I2. The rationale and code are the same as above.

Show code

# create an empty dataframe
MathStabI2 <- NULL

for (i in 1:nrow(PREDI2)) {
  
  # fetch parameter values from the equilibrium simulations df
  I1 = PREDI2$I1[i]
  I2 = PREDI2$I2[i]
  l = PREDI2$l[i]
  u1 = PREDI2$u1[i]
  u2 = PREDI2$u2[i]
  a1 = PREDI2$a1[i]
  a2 = PREDI2$a2[i]
  h1 = PREDI2$h1[i]
  h2 = PREDI2$h2[i]
  d1 = PREDI2$d1[i]
  d2 = PREDI2$d2[i]
  g = PREDI2$g[i]
  m = PREDI2$m[i]
  c = PREDI2$c[i]
  e1 = PREDI2$e1[i]
  e2 = PREDI2$e2[i]
  
  # Put in the solutions to the equilibrium values here:
  dN1 = PREDI2$N1[i]
  dP1 = PREDI2$P1[i]
  dC1 = PREDI2$C1[i]
  dN2 = PREDI2$N2[i]
  dP2 = PREDI2$P2[i]
  dC2 = PREDI2$C2[i]
  
  # create the Jacobian from Mathematica
  # NOTE: in transposing from Mathematica, the Jacobian is here assembled by row
  Jacob <- rbind(c(-l-dP1*u1, h1-dN1*u1, d1, 0, 0, 0), 
                 c(dP1*u1, -a1*dC1-h1+dN1*u1, -a1*dP1, 0, 0, 0), 
                 c(0, a1*dC1*e1, -d1-g+a1*e1*dP1, 0, 0, 0), 
                 c(0, 0, 0, -l-dP2*u2, h2-dN2*u2, d2), 
                 c(0, 0, 0, dP2*u2, -a2*dC2-h2+dN2*u2, -a2*dP2), 
                 c(0, 0, (g*m)/(c+m), 0, a2*dC2*e2, -d2+a2*e2*dP2))
  
  # Check to see if the real part of the leading eigenvalue is less than 0 or not - if it is than the system is stable.
  if (max(Re(base::eigen(Jacob)$values)) < 0 ){ 
    stable <- "stable"
  } else {
    stable <- "unstable"
  }
  
  MathStabI2 <- rbind(MathStabI2, 
                      data.frame(TIME=i, I1 , I2, l, u1, u2, a1, a2, h1, h2, 
                                 d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
                                 dP2, dC2, EiV1=Re(eigen(Jacob)$values[1]),
                                 EiV2=Re(eigen(Jacob)$values[2]),
                                 EiV3=Re(eigen(Jacob)$values[3]),
                                 EiV4=Re(eigen(Jacob)$values[4]),
                                 EiV5=Re(eigen(Jacob)$values[5]),
                                 EiV6=Re(eigen(Jacob)$values[6]),
                                 maxEv=max(Re(base::eigen(Jacob)$values)),
                                 stable=stable, biosense=PREDI2$biosense[i]))
}

MathStabI2$stable <- as.factor(MathStabI2$stable)

# separate unstable equilibria to work with later
MathStabI2US <- subset(MathStabI2, MathStabI2$stable == "unstable")

3.24% of the stability analyses runs produced unstable results (i.e., 324 out of 10000 iterations). Now we check if these are the same ones that produce equilibria where state variables are < 0 at equilibrium (i.e., not biologicaly plausible).

Show code

I2biononsense <- subset(PREDI2, PREDI2$biosense == "no")

identical(as.numeric(MathStabI2US[ , "TIME"]), I2biononsense[ , "TIME"])

[1] TRUE

Since it appears that they are, let’s add the stable/unstable information to PREDI2.

Show code

PREDI2 <- left_join(PREDI2, select(MathStabI2, !c(dN1:maxEv)), 
                    by = c("TIME", "I1", "I2", "l", "u1", "u2", "a1", "a2", "h1", 
                           "h2", "d1", "d2", "g", "m", "c", "e1", "e2", "biosense"))

Now, let’s exclude the unstable parameter sets from the analyses below and from our results. Then, we will sample 1000 iterations out of PREDI2 to use in later figures.

Show code

PREDI2pos <- filter(PREDI2, PREDI2$biosense == "yes" & PREDI2$stable == "stable")

# sample PREDpos to only use 1000 random simulation results
PREDI2_sample1000 <- PREDI2pos[sample(nrow(PREDI2pos), size = 1000), ]

PREDI2_1k <- droplevels(PREDI2_sample1000)

Graphs

Show code

# pivot the dataframe
PREDI2_biomass_long <- pivot_longer(PREDI2_1k,
                                    names_to = "Compartment",
                                    values_to = "Stock",
                                    cols = c(N1, P1, C1, N2, P2, C2, Q,
                                             Ntot, Ptot, Ctot)) %>%
  dplyr::mutate(.,
                Ecosystem = if_else(Compartment == "N1" |
                                      Compartment == "P1" |
                                      Compartment == "C1", "Donor",
                                    if_else(Compartment == "N2" |
                                              Compartment == "P2" |
                                              Compartment=="C2",
                                            "Recipient",
                                            if_else(Compartment =="Q",
                                                    "Dispersers' Pool",
                                                    "Meta-ecosystem")))) %>%
  dplyr::mutate(.,
                Ecosystem = fct_relevel(Ecosystem,
                                        "Donor",
                                        "Recipient",
                                        "Meta-ecosystem")) %>%
  dplyr::filter(.,
                Ecosystem != "Dispersers' Pool")

PREDI2_biomass_long$Compartment <- as_factor(PREDI2_biomass_long$Compartment)
PREDI2_biomass_long$Ecosystem <- as_factor(PREDI2_biomass_long$Ecosystem)

comparts_medians <- PREDI2_biomass_long %>%
  dplyr::group_by(Compartment, Ecosystem) %>%
  dplyr::summarise(., median = median(Stock), .groups = "keep")

I2bm <- ggplot(data=PREDI2_biomass_long,
               aes(x = Compartment, y = Stock, fill = Ecosystem, col = Ecosystem)) +
  geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE, outlier.shape = 25, jitter.shape = 21) +
  geom_label(data = comparts_medians,
             aes(y = median, label = round(median,2)),
             col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
             nudge_x = -0.18, label.padding = unit(0.15, "lines"),
             label.size = 0.15) +
  facet_grid(.~Ecosystem, scales = "free") +
  scale_color_manual(values = met.brewer("Egypt", 4)) +
  scale_fill_manual(values = met.brewer("Egypt", 4)) +
  ggtitle("Organic Biomass and Nutrients Stock") +
  theme_pubr() +
  coord_cartesian(y = c(0, 100)) +
  theme(legend.position = "none")

I2bm

Figure 6: Effects of along-gradient, diffusive consumers movement in a meta-ecosystem. Compared to the control scenario, the release of primary producers in ecosystem 1 is larger following movement of consumers out of this ecosystem (confront this figure with Figure 1). Environmental fertility values: I₁ = 18, I₂ = 2. All other specifications as in Figure 1.

Show code

FluxPREDI2 <- PREDI2_1k %>%
  dplyr::select(.,
                TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
  tidyr::pivot_longer(.,
                      names_to = "Compartment",
                      values_to = "Flux",
                      cols = c(FLUX_P1:MetaEcoFlux)) %>%
  dplyr::mutate(.,
                Scale = if_else(Compartment == "FLUX_P1"
                                | Compartment == "FLUX_C1"
                                | Compartment == "FLUX_Eco_1",
                                "Donor", if_else(Compartment == "FLUX_P2" |
                                                   Compartment == "FLUX_C2" |
                                                   Compartment == "FLUX_Eco_2",
                                                 "Recipient",
                                                 "Meta-ecosystem"))) %>%
  dplyr::mutate(.,
                Scale = fct_relevel(Scale,
                                       "Donor",
                                       "Recipient",
                                       "Meta-ecosystem"))

FluxPREDI2$Scale <- as.factor(FluxPREDI2$Scale)
FluxPREDI2$Compartment <- as.factor(FluxPREDI2$Compartment)

comparts_medians <- FluxPREDI2 %>%
  dplyr::group_by(Compartment, Scale) %>%
  dplyr::summarise(., median = median(Flux), .groups = "keep")

I2nf <- FluxPREDI2 %>%
  dplyr::mutate(.,
                Compartment = fct_relevel(Compartment,
                                          c("FLUX_P1", "FLUX_C1", "FLUX_Eco_1",
                                            "FLUX_P2", "FLUX_C2", "FLUX_Eco_2",
                                            "MetaEcoFlux"))) %>%
  ggplot(aes(x = Compartment, y = Flux)) +
  geom_boxjitter(outlier.intersect = TRUE, alpha = 0.1, outlier.shape = 25, jitter.shape = 21) +
  geom_label(data = comparts_medians,
             aes(y = median, label = round(median,2)),
             col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75,
             nudge_x = -0.18, label.padding = unit(0.15, "lines"),
             label.size = 0.15) +
  facet_grid(.~Scale, scales = "free") +
  coord_cartesian(y = c(0, 300)) +
  scale_x_discrete(labels = c("FLUX_C1" = "C1",
                              "FLUX_P1" = "P1",
                              "FLUX_Eco_1" = "Donor",
                              "FLUX_C2" = "C2",
                              "FLUX_P2" = "P2",
                              "FLUX_Eco_2" = "Recipient",
                              "MetaEcoFlux" = "Meta-ecosystem")) +
  ggtitle("Nutrient Flux") +
  theme_pubr() +
  theme(legend.position = "none", axis.title.x = element_blank())

I2nf

Figure 7: With a higher fertility in the ecosystem 1 (i.e., I₁ >> I₂) nutrient flux in this ecosystem increases, most likely as a consequence of the release primary producers. Conversely, ecosystem 2 shows lower and less variable levels of nutrient flux. At the meta-ecosystem scale, nutrient flux is higher than in either local ecosystems but lower than in the case of equal environmental fertility across ecosystems (Figure 5). Environmental fertility values: I₁ = 18, I₂ = 2. All other specifications as in Figure 5.

Show code

# I2all <- I2bm + I2nf + plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") + plot_layout(ncol = 1, nrow = 2)

#ggsave(I2all, filename = "../Results/LowI2BmNf.pdf", device = "pdf", dpi = 300, width = 8, height = 6)

Against-gradient consumer movement

The code chunk below analyses the experimental scenario of higher environmental fertility in ecosystem 2. In this case, consumer movement happens counter to the resource availability gradient (i.e., non-diffusive movement; Gounand et al. (2018)).

Show code

# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 < I2
I1=2
I2=18

# leaching rate
l = 0.1

# Allocate an empty data frame to save the data in
PREDI1 <- NULL
# Allocate some rows and columns. We will replace the current numbers with results below
# Data we are calculating are (bio)mass for soil nutrients (N1, N2), producers (P1, P2), consumers (C1,C2), disperses (Q)

PREDI1 <-rbind(PREDI1, 
               data.frame(TIME=seq(0,100,1), I1 = I1,  I2 = I2, l = l, 
                          u1 = seq(0,100,1), u2 = seq(0,100,1), a1 = seq(0,100,1), 
                          a2 = seq(0,100,1), h1 = seq(0,100,1), h2 = seq(0,100,1), 
                          d1 = seq(0,100,1), d2 = seq(0,100,1), g = seq(0,100,1), 
                          m = seq(0,100,1), c = seq(0,100,1), e1 = seq(0,100,1), 
                          e2 = seq(0,100,1), N1=seq(0,100,1), P1=seq(0,100,1), 
                          C1=seq(0,100,1), N2=seq(0,100,1), P2=seq(0,100,1), 
                          C2=seq(0,100,1), Q = seq(0,100,1), FLUX_P1 = seq(0,100,1), 
                          FLUX_C1 = seq(0,100,1), FLUX_P2 = seq(0,100,1), 
                          FLUX_C2 = seq(0,100,1), FLUX_Eco_1 = seq(0,100,1), 
                          FLUX_Eco_2 = seq(0,100,1), FLUX_Q = seq(0,100,1), 
                          PROD_P1 = seq(0,100,1), PROD_C1 = seq(0,100,1), 
                          PROD_P2 = seq(0,100,1), PROD_C2 = seq(0,100,1)))

# We want to start simulations at 0
i1 = 0

for (i in seq(0,9999,1)){
  # go to the next row
  i1 = i1 + 1

  # Let's sample the parameter values from DATA. We will sample these x times as defined by the i for loop above
  # producers uptake rates
  u1 = DATA1[i1,2] # in ecosystem 1
  u2 = DATA1[i1,3] # in ecosystem 2
  # consumers attack rates
  a1 = DATA1[i1,4]
  a2 = DATA1[i1,5]
  # death rates (= recycling rate)
  h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
  h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystem 2
  d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
  d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
  c = DATA1[i1,10] # death rate in the disperser's pool Q
  # movement rate to the Disperser's pool Q
  g = DATA1[i1,11]
  # movement rate from the Disperser's pool Q
  m = DATA1[i1,12]
  # efficiency of consumers 
  e1 = DATA2[i1,2] # in ecosystem 1
  e2 = DATA2[i1,3] # in ecosystem 2
  
  PREDI1[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
                 ((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (d1 + g)/(a1*e1),
                 (e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
                 (l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 ((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
                 (e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 ((d1 + g)*h1)/(a1*e1),
                 (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
                 (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 ((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 (h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
                 (c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 ((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
                 (e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
                 ((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l -I1*u1)))*u2)),
                 (e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
                 )
}

# calculate meta-ecosystem biomass stocks for each trophic compartment,
# a.k.a., meta-ecosystem productivity
PREDI1$Ntot <- PREDI1$N1 + PREDI1$N2
PREDI1$Ptot <- PREDI1$P1 + PREDI1$P2
PREDI1$Ctot <- PREDI1$C1 + PREDI1$C2

# now calculate recycling flux for the local ecosystem
PREDI1$FLUX_Eco1_check <- PREDI1$FLUX_P1 + PREDI1$FLUX_C1
PREDI1$FLUX_Eco2_check <- PREDI1$FLUX_P2 + PREDI1$FLUX_C2
# and the meta-ecosystem recycling flux
PREDI1$FLUX_Ptot <- PREDI1$FLUX_P1 + PREDI1$FLUX_P2
PREDI1$FLUX_Ctot <- PREDI1$FLUX_C1 + PREDI1$FLUX_C2
PREDI1$MetaEcoFlux <- PREDI1$FLUX_Eco_1 + PREDI1$FLUX_Eco_2 - PREDI1$FLUX_Q

# Finally, calculate the compartment production at meta-ecosystem scale
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")

# Equilibrium stocks less than or equal to 0 make no biological sense, so we 
# will remove them from the PRED dataset.

PREDI1 <- PREDI1 %>% 
  dplyr::mutate(., biosense = ifelse(N1 > 0 & 
                                     P1 > 0 & 
                                     C1 > 0 & 
                                     N2 > 0 & 
                                     P2 > 0 & 
                                     C2 > 0 & 
                                     Q > 0, 
                                    "yes", 
                                    "no"), 
                .after = MetaEcoFlux) %>% 
  dplyr::mutate_at(vars(biosense), factor)

Stability Analyses

Here we perform stability analyses for the scenario of higher environmental fertility in ecosystem 2, the recipient (i.e., I1 << I2). The analyses follow the established workflow (see above).

Show code

# create an empty dataframe
MathStabI1 <- NULL

for (i in 1:nrow(PREDI1)) {
  
  # fetch parameter values from the equilibrium simulations df
  I1 = PREDI1$I1[i]
  I2 = PREDI1$I2[i]
  l = PREDI1$l[i]
  u1 = PREDI1$u1[i]
  u2 = PREDI1$u2[i]
  a1 = PREDI1$a1[i]
  a2 = PREDI1$a2[i]
  h1 = PREDI1$h1[i]
  h2 = PREDI1$h2[i]
  d1 = PREDI1$d1[i]
  d2 = PREDI1$d2[i]
  g = PREDI1$g[i]
  m = PREDI1$m[i]
  c = PREDI1$c[i]
  e1 = PREDI1$e1[i]
  e2 = PREDI1$e2[i]
  
  # Put in the solutions to the equilibrium values here:
  dN1 = PREDI1$N1[i]
  dP1 = PREDI1$P1[i]
  dC1 = PREDI1$C1[i]
  dN2 = PREDI1$N2[i]
  dP2 = PREDI1$P2[i]
  dC2 = PREDI1$C2[i]
  
  # create the Jacobian from Mathematica
  # NOTE: in transposing from Mathematica, the Jacobian is here assembled by row
  Jacob <- rbind(c(-l-dP1*u1, h1-dN1*u1, d1, 0, 0, 0), 
                 c(dP1*u1, -a1*dC1-h1+dN1*u1, -a1*dP1, 0, 0, 0), 
                 c(0, a1*dC1*e1, -d1-g+a1*e1*dP1, 0, 0, 0), 
                 c(0, 0, 0, -l-dP2*u2, h2-dN2*u2, d2), 
                 c(0, 0, 0, dP2*u2, -a2*dC2-h2+dN2*u2, -a2*dP2), 
                 c(0, 0, (g*m)/(c+m), 0, a2*dC2*e2, -d2+a2*e2*dP2))
  
  # Check to see if the real part of the leading eigenvalue is less than 0 or not - if it is than the system is stable.
  if (max(Re(base::eigen(Jacob)$values)) < 0 ){ 
    stable <- "stable"
  } else {
    stable <- "unstable"
  }
  
  MathStabI1 <- rbind(MathStabI1, 
                      data.frame(TIME=i, I1 , I2, l, u1, u2, a1, a2, h1, h2, 
                                 d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
                                 dP2, dC2, EiV1=Re(eigen(Jacob)$values[1]),
                                 EiV2=Re(eigen(Jacob)$values[2]),
                                 EiV3=Re(eigen(Jacob)$values[3]),
                                 EiV4=Re(eigen(Jacob)$values[4]),
                                 EiV5=Re(eigen(Jacob)$values[5]),
                                 EiV6=Re(eigen(Jacob)$values[6]),
                                 maxEv=max(Re(base::eigen(Jacob)$values)),
                                 stable=stable, biosense=PREDI1$biosense[i]))
}

MathStabI1$stable <- as.factor(MathStabI1$stable)

# separate unstable equilibria to work with later
MathStabI1US <- subset(MathStabI1, MathStabI1$stable == "unstable")

2.85% of the stability analyses runs produced unstable results (i.e., 285 out of 10000 iterations). As always, let’s check whether these are also the parameter sets that produce state variable values at equilibrium that do not have biological meaning (i.e., <0).

Show code

I1biononsense <- subset(PREDI1, PREDI1$biosense == "no")

identical(as.numeric(MathStabI1US[ , "TIME"]), I1biononsense[ , "TIME"])

[1] TRUE

Now, since that proved to be the case, let’s add the stable/unstable information to PREDI1.

Show code

PREDI1 <- left_join(PREDI1, select(MathStabI1, !c(dN1:maxEv)), 
                    by = c("TIME", "I1", "I2", "l", "u1", "u2", "a1", "a2", "h1", 
                           "h2", "d1", "d2", "g", "m", "c", "e1", "e2", "biosense"))

Now, let’s exclude the unstable, biological nonsense parameter sets from the analyses below and from our results. Then, we sample 1000 random iterations from PREDI1 and store them for later use in figures and comparisons.

Show code

PREDI1pos <- filter(PREDI1, PREDI1$biosense == "yes" & PREDI1$stable == "stable")

# sample PREDpos to only use 1000 random simulation results
PREDI1_sample1000 <- PREDI1pos[sample(nrow(PREDI1pos), size = 1000), ]

PREDI1_1k <- droplevels(PREDI1_sample1000)

Graphs

Show code

# pivot the
PREDI1_biomass_long <- pivot_longer(PREDI1_1k, names_to = "Compartment", values_to = "Stock", cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot, Ctot)) %>% dplyr::mutate(., Ecosystem = if_else(Compartment == "N1"|Compartment == "C1"|Compartment == "P1", "Donor", if_else(Compartment == "N2"|Compartment == "P2"|Compartment=="C2", "Recipient", if_else(Compartment == "Q", "Dispersers' Pool", "Meta-ecosystem")))) %>% dplyr::mutate(., Ecosystem = fct_relevel(Ecosystem, "Donor", "Recipient", "Meta-ecosystem")) %>% dplyr::filter(., Ecosystem != "Dispersers' Pool")

PREDI1_biomass_long$Compartment <- as_factor(PREDI1_biomass_long$Compartment)
PREDI1_biomass_long$Ecosystem <- as_factor(PREDI1_biomass_long$Ecosystem)

comparts_medians <- PREDI1_biomass_long %>% dplyr::group_by(Compartment, Ecosystem) %>% dplyr::summarise(., median = median(Stock), .groups = "keep")

I1bm <- ggplot(data=PREDI1_biomass_long, aes(x = Compartment, y = Stock, fill = Ecosystem, col = Ecosystem)) +
  geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE, outlier.shape = 25, jitter.shape = 21) +
  geom_label(data = comparts_medians, aes(y = median, label = round(median,2)), col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75, nudge_x = -0.18, label.padding = unit(0.15, "lines"), label.size = 0.15) +
  facet_grid(.~Ecosystem, scales = "free") +
  scale_color_manual(values = met.brewer("Egypt", 4)) +
  scale_fill_manual(values = met.brewer("Egypt", 4)) +
  ggtitle("Organic Biomass and Nutrient Stock") +
  theme_pubr() +
  coord_cartesian(y = c(0, 100)) +
  theme(legend.position = "none")

I1bm

Figure 8: Movement of consumer against the environmental fertility gradient leads to a strong reduction of the nutrients stock in the donor ecosystem 1, stemming from the release of local primary producers from consumer control. In the recipient ecosystem 2, the nutrients stock increases beyond the level of enrichment seen in Figure 1, as local and immigrant consumer strongly suppress the local primary producers. All other specifications as in Figure 1. Environmental fertility values: I₁ = 2, I₂ = 18.

Following from the stronger trophic cascades, nutrient flux is also affected at both local and meta-ecosystem scales when consumer movement happens against the fertility gradient. Here, median nutrient flux in the recipient, and more fertile, ecosystem 2 is more than double the median flux in the case of diffusive consumer movement (Figure 7). Furthermore, it is also higher than in the case of equal local environmental fertility—i.e., lack of a fertility gradient (Figure 5). Likewise, median nutrient flux at the meta-ecosystem scale is higher than in either the along-gradient (Figure 7) and no-gradient (Figure 5) cases.

Show code

FluxPREDI1 <- PREDI1_1k %>% dplyr::select(TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>% tidyr::pivot_longer(names_to = "Compartment", values_to = "Flux", cols = c(FLUX_P1:MetaEcoFlux)) %>% dplyr::mutate(Scale = if_else(Compartment == "FLUX_P1" | Compartment == "FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment == "FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2", "Recipient", "Meta-ecosystem"))) %>% dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))

FluxPREDI1$Scale <- as.factor(FluxPREDI1$Scale)
FluxPREDI1$Compartment <- as.factor(FluxPREDI1$Compartment)

comparts_medians <- FluxPREDI1 %>% dplyr::group_by(Compartment, Scale) %>% dplyr::summarise(., median = median(Flux), .groups = "keep")

I1nf <- FluxPREDI1 %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1", "FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
  ggplot(., aes(x = Compartment, y = Flux)) +
  geom_boxjitter(outlier.intersect = TRUE, alpha = 0.1, outlier.shape = 25, jitter.shape = 21) +
  geom_label(data = comparts_medians, aes(y = median, label = round(median,2)), col = "black", size = 3, fill = "#FFFFFF", alpha = 0.75, nudge_x = -0.18, label.padding = unit(0.15, "lines"), label.size = 0.15) +
  facet_grid(.~Scale, scales = "free") +
  theme_pubr() +
  coord_cartesian(y = c(0, 400)) +
  scale_x_discrete(labels = c("FLUX_C1" = "C1", "FLUX_P1" = "P1", "FLUX_Eco_1" = "Donor", "FLUX_C2" = "C2", "FLUX_P2" = "P2", "FLUX_Eco_2" = "Recipient", "MetaEcoFlux" = "Meta-ecosystem")) +
  ggtitle("Nutrient Flux") +
  theme_pubr() +
  theme(legend.position = "none", axis.title.x = element_blank())

I1nf

Figure 9: Nutrient flux more than doubles in ecosystem 2, when this is the more fertile of the two ecosystem and is also the recipient of consumer movement. In turn, this leads to an overall increase in the nutrient flux at the meta-ecosystem scale, compared to scenarios where environmental fertility is equal among local ecosystem (Figure 5) or is higher in the donor ecosystem (Figure 19). Environmental fertility values: I₁ = 2, I₂ = 18.

Show code

# I1all <- I1bm + I1nf + plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") + plot_layout(ncol = 1, nrow = 2)

#ggsave(I1all, filename = "../Results/LowI1BmNf.pdf", device = "pdf", dpi = 300, width = 8, height = 6)

Summary of Stability Analyses

Here, we summarize the results of the model’s stability analyses, and collect this information to produce Table S1 shown in the Supplementary Materials. As the table shows, in all three scenario tested, only a small portion of the parameter sets used to simulate the model’s behavior at equilibrium produces unstable results. As detailed above, we exclude these unstable parameter sets from further investigations.

The table is saved in folder ../Results/ as a .tex file.

Show code

# summarise the number of stable and unstable parameter sets for each 
# stability analysis run with the numerical approximation method and 
# calculate the percentage of unstable parameter sets over the total 
# iterations of the simulations
MStabSum <- MathStab %>%
  dplyr::group_by(., stable) %>% 
  dplyr::summarise(., n = n(), .groups = "keep") %>% 
  add_column(., Fertility = "Equal fertility") %>% 
  pivot_wider(id_cols = Fertility, 
              names_from = stable, 
              values_from = n) %>% 
  dplyr::mutate(., Percent = ((unstable/nrow(MathStab))*100)) 

I2MStabSum <- MathStabI2 %>% 
  dplyr::group_by(., stable) %>% 
  dplyr::summarise(., n = n(), .groups = "keep") %>% 
  add_column(., Fertility = "Higher fertility in ecosystem 1") %>% 
  pivot_wider(id_cols = Fertility, 
              names_from = stable, 
              values_from = n) %>% 
  dplyr::mutate(., Percent = ((unstable/nrow(MathStabI2))*100)) 

I1MStabSum <- MathStabI1 %>% 
  dplyr::group_by(., stable) %>% 
  dplyr::summarise(., n = n(), .groups = "keep") %>% 
  add_column(., Fertility = "Higher fertility in ecosystem 2") %>% 
  pivot_wider(id_cols = Fertility, 
              names_from = stable, 
              values_from = n) %>% 
  dplyr::mutate(., Percent = ((unstable/nrow(MathStabI1))*100)) 

# Bind the objects generated above into a single dataframe and then produce a 
# table that shows the number of stable and unstable parameter sets for each 
# stability analysis method 

SAsummary <- bind_rows(MStabSum, I2MStabSum, I1MStabSum) %>% 
  ungroup() %>%
  gt() %>% 
  fmt_number(columns = 4, decimals = 2) %>% 
  cols_label(., Fertility = "Fertility", 
             stable = md("Stable"), 
             unstable = md("Unstable"), 
             Percent = md("\\%")) %>% 
  tab_header(title = md("**Summary of stability analyses**, showing the number of 
             Stable and Unstable equilibria, and the percentage of Unstable 
             equilibria over the total number of iterations (n = 10000). 
             Rows are grouped by environmental fertility gradient scenario.")) %>%
  tab_style(style = list(cell_text(weight = "bold")), 
            locations = cells_column_labels()) %>% 
  tab_style(style = cell_text(align = "left", size = 18), 
            locations = cells_title())

# save it
# gtsave(SAsummary, filename = "../Results/StabilityAnalyses_SummaryTable.tex")

# view it
SAsummary

Summary of stability analyses, showing the number of Stable and Unstable equilibria, and the percentage of Unstable equilibria over the total number of iterations (n = 10000). Rows are grouped by environmental fertility gradient scenario.
Fertility	Stable	Unstable	%
Equal fertility	9845	155	1.55
Higher fertility in ecosystem 1	9676	324	3.24
Higher fertility in ecosystem 2	9715	285	2.85

Changes to biomass distribution in local and meta-ecosystem

Here, we investigate whether different types of consumer movement between the two local ecosystems lead to changes in the distribution of biomass in the meta-ecosystem. Biomass distribution in local and meta-ecosystems can vary from bottom-heavy to top-heavy (McCauley et al. 2018). Bottom-heavy ecosystems present the classic biomass pyramid, with a large base of inorganic nutrients and primary producers supporting an increasingly smaller biomass of primary and secondary consumers. This bottom-heavy distribution is typical of terrestrial ecosystems—albeit exceptions exist (Hatton et al. 2015; McCauley et al. 2018). Conversely, in top-heavy ecosystems, a small base of inorganic nutrients and primary producers supports a large biomass of primary and secondary consumers. This inverted biomass pyramid, long seen as an exception, is increasingly recognized as widespread in certain systems—e.g., marine ecosystems (Sandin and Zgliczynski 2015; McCauley et al. 2018; Woodson, Schramski, and Joye 2018).

To investigate biomass distribution in our model, we use the Consumer to Resource biomass ratio, calculated as:

\(\displaystyle C{:}R = \frac{Biomass_{C^{\ast}}}{Biomass_{R^{\ast}}}\)

where the asterisk \(^{\ast}\) indicates the equilibrium biomass estimates. Values of \(C{:}R < 1\) identify ecosystems with a bottom-heavy biomass pyramid, whereas when \(C{:}R > 1\) the biomass pyramid is inverted and top-heavy (McCauley et al. 2018).

We begin by calculating the \(C{:}R\) ratio for each scenario investigated above. First, we create a new dataframe bmDistr, where we will store the required data: fertility conditions, movement rates values, efficiency rates values, state variable values at equilibrium, and ecosystem function values. We will also filter out those parameter sets that have no biological sense and are unstable as per the stability analyses above.

Show code

bmDistr <- bind_rows("Equal" = PRED, 
                     "Along" = PREDI2, 
                     "Against" = PREDI1, 
                     .id = "Fertility") %>% 
  select(., biosense, Fertility:I2, g:e2, N1:C2, Ntot:Ctot) %>% 
  filter(., biosense == "yes")

Now, we calculate \(C{:}R\). For our purposes, biomass values of primary producers in the local and meta-ecosystem will be the denominator of the ratio. That is, \(R = P^*_i\).

Show code

bmDistr <- dplyr::mutate(bmDistr, 
                         Eco1 = C1/P1, 
                         Eco2 = C2/P2, 
                         MetaEco = Ctot/Ptot 
                         )

Let’s look at the density distribution of these data. First, we exclude parameter sets that produce extreme values of \(C{:}R\)—that is, \(C{:}R > 10\).

Show code

ecoCR <- filter(bmDistr, Eco1 <= 10 & Eco2 <= 10 & MetaEco <= 10)
# summary(ecoCR)

Now, we pivot the resulting dataset and plot the density distribution of \(C{:}R\) values.

Show code

fert.labs <- c("Against-gradient", "Along-gradient", "Gradient-neutral")
names(fert.labs) <- c("Against", "Along", "Equal")

scale.labs <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")
names(scale.labs) <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")

ecoCR_long <- pivot_longer(ecoCR, cols = c(Eco1, Eco2, MetaEco), 
                           values_to = "CRratio", names_to = "Scale") %>% 
  mutate_at(., vars(Scale), factor) %>% 
  mutate(., Scale = ifelse(Scale == "Eco1", 
                           "Ecosystem 1", ifelse(Scale == "Eco2", 
                                                 "Ecosystem 2", 
                                                 "Meta-ecosystem")))

ggplot(ecoCR_long, aes(x = CRratio, fill = Fertility, col = Fertility)) + 
  geom_density(aes(y = ..density..), alpha = 0.1) + 
  scale_color_highcontrast(reverse = T, name = "Consumer movement", labels = c("Against-gradient", "Along-gradient", "Gradient-neutral")) +
  scale_fill_highcontrast(reverse = T, name = "Consumer movement", labels = c("Against-gradient", "Along-gradient", "Gradient-neutral")) +
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) + 
  theme_pubr(legend = "right", base_size = 10) + 
  xlab("C:R ratio") 

Figure 10: Density distribution when dataset is limited to C:R < 10. Note the different scales on the y-axis.

Show code

ggsave(filename = "../Results/Density_CRratio_lessThan10.pdf", device = "pdf", dpi = 300, width = 10, height = 7)

Movement from ecosystem 1

Here we look at how the \(C{:}R\) values change with increasing movement from ecosystem 1, the donor, towards the dispersers’ pool Q (i.e., parameter g in the model). The graph below uses the full dataset created above, ecoCR_long, and fits linear models to each combination of Fertility scenario (gradient-neutral, along-, and against-gradient) and Scale (ecosystem 1, ecosystem 2, and meta-ecosystem.

Show code

ecoCRg_lm <- ecoCR_long %>% nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ g, data = data)))  

ecoCRg_lmOutput <- ecoCRg_lm %>% summarise(glance(model))

ecoCRg_lmOutput$g_slope <- NA
ecoCRg_lmOutput$g_se <- NA

for (i in 1:length(ecoCRg_lm[["model"]])) {
  ecoCRg_lmOutput$g_slope[i] <- ecoCRg_lm$model[[i]]$coefficient[2]
  ecoCRg_lmOutput$g_se[i] <- summary(ecoCRg_lm$model[[i]])$coefficient["g", "Std. Error"]
}

ecoCRg_lmOutput %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n",
             g_slope = "g Coeff.",
             g_se = "g SE") %>% 
  gtsave(., filename = "../Results/CR_g_lm_res.rtf")

ecoCR_long %>% ggplot(., aes(x = g, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) + 
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) + 
  ylab("C:R ratio") + 
  xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values [0,10].") +
  theme_pubr(legend = "right", base_size = 10)

Figure 11: Distribution of C:R values lower than 10 as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Darker shades of blue correspond to a higher count of C:R value in a given bin. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data. C:R values appear to follow a decreasing trend as g increases in Ecosystem 1 and in the Meta-ecosystem, whereas they appear stable or slightly increasing in Ecosystem 2.

Show code

ggsave(filename = "../Results/CRratioVg_binned.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

Now, let’s explore the relationships in a bit more detail. We do this by splitting the dataframe in two portions: \(C{:}R > 1\) and \(C{:}R < 1\). Here, we take a look at the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).

Show code

ecoCRg_lm_above1 <- ecoCR_long %>% 
  filter(., CRratio > 1) %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ g, data = data)))  

ecoCRg_lmOutput_above1 <- ecoCRg_lm_above1 %>% summarise(glance(model)) 

ecoCRg_lmOutput_above1$g_slope <- NA
ecoCRg_lmOutput_above1$g_se <- NA

for (i in 1:length(ecoCRg_lm_above1[["model"]])) {
  ecoCRg_lmOutput_above1$g_slope[i] <- ecoCRg_lm_above1$model[[i]]$coefficient[2]
  ecoCRg_lmOutput_above1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g", "Std. Error"]
}

ecoCRg_lmOutput_above1 %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n", 
             g_slope = "g Coeff.",
             g_se = "g SE") %>% 
  tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers' 
                        Pool (g)"), 
             subtitle = md("C:R values between (1,10].")) %>% 
  gtsave(., filename = "../Results/CR_g_lm_above1.rtf")

ecoCR_long %>% 
  filter(., CRratio > 1) %>% 
  ggplot(., aes(x = m, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) + 
  ylab("C:R ratio") + 
  xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values between (1,10].") +
  theme_pubr(legend = "right", base_size = 10)

Figure 12: Distribution of C:R values comprised between (1, 10] as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. As g increases, we do not observe any trend in the values of C:R. Note that C:R > 1 are generally indicative of top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVg_binned_above1.pdf", device = "pdf", 
       dpi = 300, width = 12, height = 12)

And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).

Show code

ecoCRg_lm_below1 <- ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ g, data = data)))  

ecoCRg_lmOutput_below1 <- ecoCRg_lm_below1 %>% summarise(glance(model))

ecoCRg_lmOutput_below1$g_slope <- NA
ecoCRg_lmOutput_below1$g_se <- NA

for (i in 1:length(ecoCRg_lm_below1[["model"]])) {
  ecoCRg_lmOutput_below1$g_slope[i] <- ecoCRg_lm_below1$model[[i]]$coefficient[2]
  ecoCRg_lmOutput_below1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g", "Std. Error"]
}

ecoCRg_lmOutput_below1 %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(., columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n", 
             g_slope = "g Coeff.",
             g_se = "g SE") %>% 
  tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers' 
                        Pool (g)"), 
             subtitle = md("C:R values between [0,1].")) %>% 
  gtsave(., filename = "../Results/CR_g_lm_below1.rtf")

ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  ggplot(., aes(x = m, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) +  
 facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +  
  ylab("C:R ratio") + 
  xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values between [0,1].") + 
  scale_color_manual(name = "Legend", values = c("cyan", "magenta")) + 
  theme_pubr(legend = "right", base_size = 10)

Figure 13: Distribution of C:R values [0,1] when movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Again, no trend is apparent save for a seeming increase in C:R values in ecosystem 2 when consumers move along-gradient (central panel). Note, however, that this trend does not scale up to the meta-ecosystem. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVg_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

Movement towards ecosystem 2

Here, we investigate how the \(C{:}R\) biomass ratio changes with increasing values of m, the movement rate of consumers leaving the dispersers’ pool Q to enter ecosystem 2. Here is the general relationship.

Show code

ecoCRm_lm <- ecoCR_long %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ m, data = data)))  

ecoCRm_lmOutput <- ecoCRm_lm %>% summarise(glance(model))

ecoCRm_lmOutput$m_slope <- NA
ecoCRm_lmOutput$m_se <- NA

for (i in 1:length(ecoCRm_lm[["model"]])) {
  ecoCRm_lmOutput$m_slope[i] <- ecoCRm_lm$model[[i]]$coefficient[2]
  ecoCRm_lmOutput$m_se[i] <- summary(ecoCRm_lm$model[[i]])$coefficient["m", "Std. Error"]
}

ecoCRm_lmOutput %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value,logLik:deviance, m_slope:m_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n",
             m_slope = "m Coeff.",
             m_se = "m SE") %>% 
  gtsave(., filename = "../Results/CR_m_lm_res.rtf")

ecoCR_long %>% 
  ggplot(., aes(x = m, y = CRratio, )) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
 geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) + 
  ylab("C:R ratio") + 
  xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values limited to [0,10].") + 
  scale_color_manual(name = "Legend", values = c("cyan", "magenta")) +
  theme_pubr(legend = "right", base_size = 10)

Figure 14: Density distribution of C:R values [0,10] as consumer movement rate from the dispersers’ pool to ecosystem 2 (m) increases. No trend is clearly visible, aside from an apparent increase in C:R in ecosystem 2 when consumers move along-gradients (central panel). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVm_binned.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).

Show code

ecoCRm_lm_above1 <- ecoCR_long %>% 
  filter(., CRratio > 1) %>% 
  nest_by(Fertility, 
          Scale) %>% 
  mutate(model = list(lm(CRratio ~ m, 
                         data = data)))  

ecoCRm_lmOutput_above1 <- ecoCRm_lm_above1 %>% summarise(glance(model)) 

ecoCRm_lmOutput_above1$m_slope <- NA
ecoCRm_lmOutput_above1$m_se <- NA

for (i in 1:length(ecoCRm_lm_above1[["model"]])) {
  ecoCRm_lmOutput_above1$m_slope[i] <- ecoCRm_lm_above1$model[[i]]$coefficients[2]
  ecoCRm_lmOutput_above1$m_se[i] <- summary(ecoCRm_lm_above1$model[[i]])$coefficient["m", "Std. Error"]
}

ecoCRm_lmOutput_above1 %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value, 
                         logLik:deviance, 
                         m_slope:m_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n", 
             m_slope = "m Coeff.",
             m_se = "m SE") %>% 
  tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers' 
                        Pool (m)"), 
             subtitle = md("C:R values between (1,10].")) %>% 
  gtsave(., filename = "../Results/CR_m_lm_above1.rtf")

ecoCR_long %>% 
  filter(., CRratio > 1) %>% 
  ggplot(., aes(x = m, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +  
  ylab("C:R ratio") + 
  xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values between (1,10].") + 
  scale_color_manual(name = "Legend", values = c("cyan", "magenta")) +
  theme_pubr(legend = "right", base_size = 10)

Figure 15: Distribution of C:R (1, 10] as movement rate from the dispersers’ pool to ecosystem 2 (m) increases. Focusing on only C:R values indicative of inverted, top-heavy biomass pyramids shows weakly increasing trends in the value of C:R in ecosystem 1 when consumers move against-gradient (top-left panel), and in ecosystem 2 when consumers move along-gradient (central panel). We also note an apparent decrese in C:R values in ecosystem 1 as consumers move in an homogeneous meta-ecosystem (bottom-left panel). Note, however, that none of these trends appears to influence biomass dynamics at the meta-ecosystem scale (right-most column). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVm_binned_above1.pdf", device = "pdf", 
       dpi = 300, width = 12, height = 12)

And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).

Show code

ecoCRm_lm_below1 <- ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ m, data = data)))  

ecoCRm_lmOutput_below1 <- ecoCRm_lm_below1 %>% summarise(glance(model))

ecoCRm_lmOutput_below1$m_slope <- NA
ecoCRm_lmOutput_below1$m_se <- NA

for (i in 1:length(ecoCRm_lm_below1[["model"]])) {
  ecoCRm_lmOutput_below1$m_slope[i] <- ecoCRm_lm_below1$model[[i]]$coefficient[2]
  ecoCRm_lmOutput_below1$m_se[i] <- summary(ecoCRm_lm_below1$model[[i]])$coefficient["m", "Std. Error"]
}

ecoCRm_lmOutput_below1 %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value,logLik:deviance, m_slope:m_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n", 
             m_slope = "m Coeff.",
             m_se = "m SE") %>% 
  tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers' 
                        Pool (m)"), 
             subtitle = md("C:R values between [0,1].")) %>% 
  gtsave(., filename = "../Results/CR_m_lm_below1.rtf")

ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  ggplot(., aes(x = m, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) + 
  ylab("C:R ratio") + 
  xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values between [0,1].") + 
  theme_pubr(legend = "right", base_size = 10)

Figure 16: Distribution of C:R [0,1] as consumer movement rate from the dispersers pool to ecosystem 2 (m) increases. Here, we notice apparently increasing trends in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). We do not see any influence of these trends at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVm_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

Consumer death rate while in transit

Here, we investigate how the \(C{:}R\) biomass ratio changes with increasing values of c, the death rate of consumers while in the dispersers’ pool. Consumer biomass lost while transiting through Q is lost from the meta-ecosystem and does not contribute to the recycling pathways at either local or regional scale. Hence, it is important to track how biomass distribution in the system may vary with this parameter. Here is the general relationship.

Show code

ecoCRc_lm <- ecoCR_long %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ c, data = data)))  

ecoCRc_lmOutput <- ecoCRc_lm %>% summarise(glance(model))

ecoCRc_lmOutput$c_slope <- NA
ecoCRc_lmOutput$c_se <- NA

for (i in 1:length(ecoCRc_lm[["model"]])) {
  ecoCRc_lmOutput$c_slope[i] <- ecoCRc_lm$model[[i]]$coefficient[2]
  ecoCRc_lmOutput$c_se[i] <- summary(ecoCRc_lm$model[[i]])$coefficient["c", "Std. Error"]
}

ecoCRc_lmOutput %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value,logLik:deviance, c_slope:c_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n",
             c_slope = "c Coeff.",
             c_se = "c SE") %>% 
  gtsave(., filename = "../Results/CR_c_lm_res.rtf")

ecoCR_long %>% 
  ggplot(., aes(x = c, y = CRratio, )) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +  
  ylab("C:R ratio") + 
  xlab("Consumers death rate while in the Dispersers' Pool (c)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values limited to [0,10].") + 
  theme_pubr(legend = "right", base_size = 10)

Figure 17: Distribution of C:R values [0,10] as the consumer death rate in the dispersers’ pool (c) increases. No trend is apparent, aside from a weak decrease in ecosystem 2 (central and center-bottom panels) that have no influence on meta-ecosystem biomass dynamics. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVc_binned.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).

Show code

ecoCRc_lm_above1 <- ecoCR_long %>% 
  filter(., CRratio > 1) %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ c, data = data)))  

ecoCRc_lmOutput_above1 <- ecoCRc_lm_above1 %>% summarise(glance(model)) 

ecoCRc_lmOutput_above1$c_slope <- NA
ecoCRc_lmOutput_above1$c_se <- NA

for (i in 1:length(ecoCRc_lm_above1[["model"]])) {
  ecoCRc_lmOutput_above1$c_slope[i] <- ecoCRc_lm_above1$model[[i]]$coefficients[2]
  ecoCRc_lmOutput_above1$c_se[i] <- summary(ecoCRc_lm_above1$model[[i]])$coefficient["c", "Std. Error"]
}

ecoCRc_lmOutput_above1 %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value,logLik:deviance, c_slope:c_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n", 
             c_slope = "c Coeff.",
             c_se = "c SE") %>% 
  tab_header(title = md("Modeling C:R ratio vs death rate while in the
                        Dispersers' Pool (c)"), 
             subtitle = md("C:R values between (1,10].")) %>% 
  gtsave(., filename = "../Results/CR_c_lm_above1.rtf")

ecoCR_long %>% 
  filter(., CRratio > 1) %>% 
  ggplot(., aes(x = c, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) + 
  ylab("C:R ratio") + 
  xlab("Consumer death rate while in the Dispersers' Pool (c)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values between (1,10].") + 
  theme_pubr(legend = "right", base_size = 10)

Figure 18: Distribution of C:R (1,10] as consumer death rate in the dispersers’ pool (c) increases. Regardless of the type of consumer movement, we see no trend in the distribution of C:R values that identify top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVc_binned_above1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).

Show code

ecoCRc_lm_below1 <- ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  nest_by(Fertility, Scale) %>% 
  mutate(model = list(lm(CRratio ~ c, data = data)))  

ecoCRc_lmOutput_below1 <- ecoCRc_lm_below1 %>% summarise(glance(model))

ecoCRc_lmOutput_below1$c_slope <- NA
ecoCRc_lmOutput_below1$c_se <- NA

for (i in 1:length(ecoCRc_lm_below1[["model"]])) {
  ecoCRc_lmOutput_below1$c_slope[i] <- ecoCRc_lm_below1$model[[i]]$coefficient[2]
  ecoCRc_lmOutput_below1$c_se[i] <- summary(ecoCRc_lm_below1$model[[i]])$coefficient["c", "Std. Error"]
}

ecoCRc_lmOutput_below1 %>% 
  gt(groupname_col = "Scale") %>% 
  fmt_number(columns = c(r.squared:p.value,logLik:deviance, c_slope:c_se), 
             decimals = 3) %>% 
  cols_label(., r.squared = "R^2", 
             adj.r.squared = "Adj. R^2", 
             sigma = "SE", 
             df.residual = "Residual df", 
             nobs = "n", 
             c_slope = "c Coeff.",
             c_se = "c SE") %>% 
  tab_header(title = md("Modeling C:R ratio vs death rate while in the 
                        Dispersers' pool (c)"), 
             subtitle = md("C:R values between [0,1].")) %>% 
  gtsave(., filename = "../Results/CR_c_lm_below1.rtf")

ecoCR_long %>% 
  filter(., CRratio <= 1) %>% 
  ggplot(., aes(x = c, y = CRratio)) + 
  geom_bin2d(bins = 70) + 
  scale_fill_met_c("Hokusai2", direction = 1) +
  geom_smooth(method = "lm", se = FALSE, col = "black",
              linetype = 2, show.legend = T, size = .5) + 
  geom_quantile(method = "rq", quantiles = .5, col = "black",
                show.legend = T) + 
  facet_grid(Fertility~Scale, scales = "free", labeller = labeller(Fertility = fert.labs, Scale = scale.labs)) +  
  ylab("C:R ratio") + 
  xlab("Consumer death rate while in the Dispersers' Pool (c)") + 
  # ggtitle("Binned distribution of C:R", 
  #         subtitle = "C:R values between [0,1].") + 
  theme_pubr(legend = "right", base_size = 10)

Figure 19: Distribution of C:R [0,1] as consumer death rate in the dispersers’ pool (c) increases. Focusing on C:R < 1 brings out a weak, negative trend in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). Neither of these trends have any influence at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.

Show code

ggsave(filename = "../Results/CRratioVc_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)

Visual Deliverables

Here we summarize the above results and produce Figures 2, 3, and 4 shown in the main text and Figures S1 to S5 in the Supplementary Materials.

We present figures for both the raw data and the log response ratio (\(LRR_X\)) of nutrient stock and biomass accumulation, nutrient flux, and trophic compartment productivity. Our reference for calculating the \(LRR_X\) will be the results for the control scenario (Control scenario: gradient-neutral consumer movement), where environmental fertility is homogeneous between local ecosystems (i.e., I₁ = I₂ = 10).

Mathematically, the ratio’s formula is:

\(\displaystyle LRR_X = log_{10} \biggl(\frac{X_{i, I_1 \gtrless I_2}}{X_{i, I_1 = I_2}}\biggr)\)

where \(X \in [N, P, C]\) is the trophic compartment of interest, \(i\) is either the local or meta-ecosystem, and \(I_1 \gtrless I_2\) represents the direction of fertility gradient in the scenario of interest.

A value of \(LRR_X = 0\) indicates that the ecosystem function values in the experimental scenario does not differ from that in the control scenario. Values above 0 (white background) indicate that the ecosystem function values in the experimental scenario are larger than in the control—hence, consumer movement increases them. Viceversa, values of \(LRR_X < 0\) (grey background) indicate that ecosystem function values in the control scenario are larger than in the experimental one—hence, consumer movement decreases them.

All stable parameter sets

Show code

# Here we create a custom ggplot2 theme to use in the graphs below
consmov_theme <- theme_pubr() + theme(legend.position = "bottom", 
                                      legend.text = element_text(size = 16), 
                                      legend.title = element_text(size = 17), 
                                      plot.title = element_text(size = 17),
                                      axis.title.x = element_text(size = 17, 
                                                                  face = "bold"),
                                      axis.text.x = element_text(size = 15),
                                      axis.text.y = element_text(size = 15),
                                      axis.title.y = element_text(size = 17, 
                                                                  face = "bold"),
                                      strip.text = element_text(size = 16))

We are working with the PRED, PREDI2, and PREDI1 datasets we produced above (sections Control scenario: gradient-neutral consumer movement and Experimental scenarios: along- and against-gradient consumer movement).

First, we will check whether there is overlap among the equilibria that are either unstable or without biological meaning across PRED, PREDI2, and PREDI1. Then, we will join them into a new object, EnvFert_res.

Show code

PREDusns <- filter(PRED, biosense == "no" | stable == "unstable")

PREDI2usns <- filter(PREDI2, biosense == "no" | stable == "unstable")

PREDI1usns <- filter(PREDI1, biosense == "no" | stable == "unstable")

summary(PREDusns[["TIME"]] %in% PREDI1usns[["TIME"]])

   Mode   FALSE    TRUE 
logical      61      94 

Show code

summary(PREDusns[["TIME"]] %in% PREDI2usns[["TIME"]])

   Mode   FALSE    TRUE 
logical      29     126 

Show code

EnvFert_res <- bind_rows("Equal" = PRED, "Along" = PREDI2, "Against" = PREDI1, .id = "Fertility") %>% select(., biosense, Fertility:I2, N1:C2, Ntot:Ctot, FLUX_P1:FLUX_C2, FLUX_Ptot:FLUX_Ctot, PROD_P1:PROD_Ctot, u1:e2)

There is not overlap among the three datasets in terms of equilibria that are unstable and lack biological meaning. This is something that will need accounting for later on when we work with \(LRR_X\).

Note that EnvFert_respos has the same number of rows as bmDistr created in Changes to biomass distribution in local and meta-ecosystem.

For now, we can remove all equilibria that are either unstable or do not have biological sense as we are going to create a figure using the raw data.

Show code

EnvFert_respos <- filter(EnvFert_res, biosense == "yes")

Untransformed results

We now transform EnvFert_respos to a long-format dataframe and add two categorical variables: Function and Scale. Function tells us which ecosystem function is measured. Scale shows instead whether the function is being measured at the local or meta-ecosystem scale.

Show code

EnvFert_long <- select(EnvFert_respos, !u1:e2) %>% 
  pivot_longer(., cols = N1:PROD_Ctot,
               names_to = "EcoFunction", 
               values_to = "Value") %>% 
  dplyr::mutate(., Function = if_else(EcoFunction == "N1" | 
                                        EcoFunction == "N2" | 
                                        EcoFunction == "P1" | 
                                        EcoFunction == "P2" | 
                                        EcoFunction == "C1" | 
                                        EcoFunction == "C2" | 
                                        EcoFunction == "Ntot" | 
                                        EcoFunction == "Ptot" | 
                                        EcoFunction == "Ctot", 
                                      "Stock", 
                                      if_else(EcoFunction == "FLUX_P1" | 
                                                EcoFunction == "FLUX_P2" | 
                                                EcoFunction == "FLUX_C1" | 
                                                EcoFunction == "FLUX_C2" | 
                                                EcoFunction == "FLUX_Ptot" | 
                                                EcoFunction == "FLUX_Ctot", 
                                              "Nutrient Flux", 
                                              "Trophic Productivity"))) %>% 
  dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | 
                                     EcoFunction == "P1" | 
                                     EcoFunction == "C1" | 
                                     EcoFunction == "FLUX_P1" | 
                                     EcoFunction == "FLUX_C1" | 
                                     EcoFunction == "PROD_P1" | 
                                     EcoFunction == "PROD_C1", 
                                   "Ecosystem 1", 
                                   if_else(EcoFunction == "N2" | 
                                             EcoFunction == "P2" | 
                                             EcoFunction == "C2" | 
                                             EcoFunction == "FLUX_P2" | 
                                             EcoFunction == "FLUX_C2" | 
                                             EcoFunction == "PROD_P2" | 
                                             EcoFunction == "PROD_C2", 
                                           "Ecosystem 2", 
                                           "Meta-ecosystem"))) %>% 
  dplyr::mutate(., Fertility = factor(Fertility), 
                Function = factor(Function), 
                Scale = factor(Scale), 
                EcoFunction = factor(EcoFunction))

head(EnvFert_long)

# A tibble: 6 × 9
  biosense Fertility  TIME    I1    I2 EcoFunction Value Function
  <fct>    <fct>     <dbl> <dbl> <dbl> <fct>       <dbl> <fct>   
1 yes      Equal         1    10    10 N1           3.07 Stock   
2 yes      Equal         1    10    10 P1           1.43 Stock   
3 yes      Equal         1    10    10 C1           1.31 Stock   
4 yes      Equal         1    10    10 N2           3.26 Stock   
5 yes      Equal         1    10    10 P2           1.55 Stock   
6 yes      Equal         1    10    10 C2           2.68 Stock   
# … with 1 more variable: Scale <fct>

Now that we have a dataset ready for plotting, we produce the graphs. The next code chunk produces Figure S1 in the Supplementary Materials. Note that we show the code, but not the figure itself as it is presented in the Supplementary Materials.

Show code

StockSumm <- EnvFert_long %>% filter(., Function == "Stock") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1", Nutrients = "N2", Nutrients = "Ntot", "Primary\n Producers" = "P1", "Primary\n Producers" = "P2", "Primary\n Producers" = "Ptot", Consumers = "C1", Consumers = "C2", Consumers = "Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients", "Primary\n Producers", "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 10, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) +
  scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral", "Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Stock and Biomass") +
  ylab(" ") +
  xlab(" ") +
  coord_cartesian(ylim = c(0, 75)) 

FluxSumm <- EnvFert_long %>% filter(., Function == "Nutrient Flux") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, "Primary\n Producers" = "FLUX_P1", "Primary\n Producers" = "FLUX_P2", "Primary\n Producers" = "FLUX_Ptot", Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers", "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) + 
  scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral", "Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Flux") +
  ylab(" ") +
  xlab(" ") +
  coord_cartesian(ylim = c(0, 350))

ProdSumm <- EnvFert_long %>% 
  filter(., Function == "Trophic Productivity") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, "Primary\n Producers" = "PROD_P1", "Primary\n Producers" = "PROD_P2", "Primary\n Producers" = "PROD_Ptot", Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers", "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) + 
 scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral", "Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Trophic Compartment Productivity") +
  ylab(" ") +
  xlab("Trophic compartment") +
  coord_cartesian(ylim = c(0, 500)) 

FuncAll <- StockSumm + FluxSumm + ProdSumm + 
  plot_layout(ncol = 1, nrow = 3, guides = "collect") + 
  plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") & 
  consmov_theme

FuncAll

Show code

ggsave(FuncAll, filename = "../Results/FunctionSummary_10kN.pdf", device = "pdf", dpi = 300, width = 11, height = 10)

Log-transformed Response Ratio

Now, we produce Figures 2, 3, and 4 from the main text. First, we are going to create new dataframes (PREDI2rr and PREDI1rr1) that contain only the response ratio values for nutrient stock and biomass accumulation, nutrient flux, and trophic compartment productivity of all state variables. These dataframes will also contain information about the fertility conditions under which we measure each ecosystem function. After calculating the response ratio, we will remove the unstable, biologically not meaningful equilibria from PREDI1 and PREDI2 before merging them in a new object, EnvFertRR.

In calculating \(LRR_X\), the discrepancy in the number of equilibria that are either unstable or lack biological meaning among PRED, PREDI1, and PREDI2 matters. To account for this, let us first create two objects, I1NsUsrm and I2NsUsrm that contain the unstable, nonsensical equilibria that are found in PRED but are not included in either PREDI1 and PREDI2. We are going to remove the ones in PREDI1 and PREDI2 making use of the columns biosense and stable. Conversely, we will use I1NsUsrm and I2NsUsrm to remove instance of the response ratio being calculated using a denominator that is from an unstable, nonsensical equilibrium.

Show code

I1NsUsrm <- subset(PREDusns, !(PREDusns$TIME %in% PREDI1usns$TIME))

I2NsUsrm <- subset(PREDusns, !(PREDusns$TIME %in% PREDI2usns$TIME))

Then, let’s create two new objects, PREDI2rr and PREDI1rr, to store the response ratio values calculated by dividing ecosystem function values in PREDI2 and PREDI1 by those in PRED, respectively. We will then merge these two objects in a single one called EnvFertRR, and calculated the \(C{:}R\) biomass ratio.

Show code

# need to keep raw data in the same df as the ratio to select the >0 cases 
# for state variables

PREDI2rr <- PREDI2 %>%  
  # we are going to batch-create new response ratio versions of each variable
  # we start with the "stock" group of variables
  dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1, N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot, Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
  # next, we work on the "flux" group of variables
  dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1, FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2, FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
  # finally, we work on the "productivity" group of variables
  dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1, PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2, PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
  # and now let's subset the dataframe to include the resp ratio columns
  dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
  dplyr::filter(., biosense == "yes" | stable == "yes") %>%
  dplyr::filter(., !(TIME %in% I2NsUsrm$TIME))

# repeat, but for PREDI1

PREDI1rr <- PREDI1 %>%  
  # we are going to batch-create new response ratio versions of each variable
  # we start with the "stock" group of variables
  dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1, N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot, Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
  # next, we work on the "flux" group of variables
  dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1, FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2, FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
  # finally, we work on the "productivity" group of variables
  dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1, PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2, PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
  # and now let's subset the dataframe to include only the resp ratio columns
  dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
  dplyr::filter(., biosense == "yes" | stable == "yes") %>%
  dplyr::filter(., !(TIME %in% I1NsUsrm$TIME))

EnvFertRR <- bind_rows("Along" = PREDI2rr, 
                       "Against" = PREDI1rr, 
                       .id = "Fertility") %>% 
  dplyr::select(., Fertility, 
                biosense, 
                stable, 
                TIME:I2, 
                N1rr:PROD_Crr, 
                N1:Ctot, u1:e2) %>% 
  dplyr::mutate(., 
                CR_Eco1 = C1/P1, 
                CR_Eco2 = C2/P2, 
                CR_MetaEco = Ctot/Ptot 
                ) 
# now build the graph, first by pivoting to longer format, then graphing

Note that object EnvFertRR contains 19301 and not 20 000 as would be expected from merging PREDI1 and PREDI2, because of the two rounds of removal of unstable equilibria after said merge. The first round removes those unstable, biologically not meaningful equilibria contained in PREDI1 and PREDI2, the second uses objects I1NsUsrm and I2NsUsrm to remove the ones that are contained in PRED but do not overlap with those in PREDI1 or PREDI2. T

Now, we build the graphs. First, pivot the EnvFertRR dataframe from wide to long format.

Show code

EnvFertRR_long <- select(EnvFertRR, !(N1:e2)) %>% 
  pivot_longer(., 
               cols = N1rr:PROD_Crr, 
               names_to = "Compartment",
               values_to = "Ratio") %>% 
  dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | 
                                        Compartment == "FLUX_P2rr" | 
                                        Compartment == "FLUX_C1rr" | 
                                        Compartment == "FLUX_C2rr" | 
                                        Compartment == "FLUX_Prr" | 
                                        Compartment == "FLUX_Crr", 
                                      "Nutrient Flux", 
                                      if_else(Compartment == "PROD_P1rr" | 
                                                Compartment == "PROD_C1rr" | 
                                                Compartment == "PROD_P2rr" | 
                                                Compartment == "PROD_C2rr" | 
                                                Compartment == "PROD_Prr" | 
                                                Compartment == "PROD_Crr", 
                                              "Trophic Productivity", 
                                              "Stock"))) %>% 
  dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | 
                                     Compartment == "P1rr" | 
                                     Compartment == "C1rr" | 
                                     Compartment == "FLUX_P1rr" | 
                                     Compartment == "FLUX_C1rr" | 
                                     Compartment == "PROD_P1rr" | 
                                     Compartment == "PROD_C1rr", 
                                   "Ecosystem 1", 
                                   if_else(Compartment == "N2rr" | 
                                             Compartment == "P2rr" | 
                                             Compartment == "C2rr" | 
                                             Compartment == "FLUX_P2rr" | 
                                             Compartment == "FLUX_C2rr" | 
                                             Compartment == "PROD_P2rr" | 
                                             Compartment == "PROD_C2rr", 
                                           "Ecosystem 2", 
                                           "Meta-ecosystem"))) %>% 
  dplyr::mutate(., Fertility = factor(Fertility), 
                Function = factor(Function), 
                Scale = factor(Scale), 
                Compartment = factor(Compartment)) %>%
  dplyr::mutate(., biosense = fct_drop(biosense), 
                stable = fct_drop(stable))

In the following code chunks, we subset and then plot each ecosystem function individually against the LRR to create Figure 2, 3, and 4 in the main text. Note that we show the code here but not the figures themselves, referring the reader to the main text to view the resulting figures.

This code chunk produces Figure 2, showing changes in nutrient stocks and biomass accumulation following consumer movement.

Show code

ptHCsubset <- c("#DDAA33", "#BB5566")
ptHCsubset_fill <- c("#f8eed6", "#f1dde0")

EFstockLRR <- EnvFertRR_long %>% filter(., Function == "Stock") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr", Nutrients = "N2rr", Nutrients = "Nrr", "Primary\n Producers" = "P1rr", "Primary\n Producers" = "P2rr", "Primary\n Producers" = "Prr", Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients", "Primary\n Producers", "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Stock and Biomass") +
  ylab("LRR") +
  xlab("Trophic compartment") +
  consmov_theme

ggsave(EFstockLRR, filename = "../Results/Stock_LRR_10kN.pdf", device = "pdf", dpi = 300, width = 10, height = 6)

The following code chunk produces Figure 3, which shows changes to nutrient flux following consumer movement.

Show code

EFfluxLRR <- EnvFertRR_long %>% filter(., Function == "Nutrient Flux") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Primary\n Producers" = "FLUX_P1rr", "Primary\n Producers" = "FLUX_P2rr", "Primary\n Producers" = "FLUX_Prr", Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers", "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Flux") +
  ylab("LRR") +
  xlab("Trophic compartment") +
  consmov_theme

ggsave(EFfluxLRR, filename = "../Results/Flux_LRR_10kN.pdf", device = "pdf", dpi = 300, width = 10, height = 6)

And, finally, here we produce Figure 4, that shows change in the trophic compartment productivity following consumer movement.

Show code

EFprodLRR <- EnvFertRR_long %>% 
  filter(., Function == "Trophic Productivity") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Primary\n Producers" = "PROD_P1rr", "Primary\n Producers" = "PROD_P2rr", "Primary\n Producers" = "PROD_Prr", Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers", "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Trophic Compartment Productivity") +
  ylab("LRR") +
  xlab("Trophic compartment") + 
  consmov_theme

ggsave(EFprodLRR, filename = "../Results/Production_LRR_10kN.pdf", device = "pdf", dpi = 300, width = 10, height = 6)

In the two sections below, we investigate if our results are robust to changes in biomass distribution within the system caused by moving consumers and produce Figures S2 to S5 in the Supplementary Materials.

Consumer movement in bottom-heavy systems

Here, we focus on those parameter sets that produce \(C{:}R < 1\). That is, on those parameter sets that produce classic biomass pyramids (McCauley et al. 2018).

As before, we begin with graphs for the untransformed data and then proceed to the \(LRR_X\) graphs.

Untransformed parameter sets

First, using the object EnvFert_CRbelow10 created above that exclude all \(C{:}R\) values above 10, we further exclude all \(C{:}R\) values above 1 and create a new object, EnvFert_CRbelow1.

Show code

EnvFert_CRbelow1 <- dplyr::mutate(EnvFert_respos,
                         CR_Eco1 = C1/P1, # C:R for Ecosystem 1
                         CR_Eco2 = C2/P2, # C:R for Ecosystem 2
                         CR_MetaEco = Ctot/Ptot # C:R for meta-ecosystem
                         ) %>%
  filter(., CR_Eco1 <= 1 & CR_Eco2 <= 1 & CR_MetaEco <= 1)

Now, we lengthen this new dataset for plotting with ggplot2.

Show code

EnvFert_CRbelow1_long <- select(EnvFert_CRbelow1, !u1:e2) %>%
  pivot_longer(., cols = N1:CR_MetaEco, 
               names_to = "EcoFunction", 
               values_to = "Value") %>% 
  dplyr::mutate(., Function = if_else(EcoFunction == "N1" | 
                                        EcoFunction == "N2" | 
                                        EcoFunction == "P1" | 
                                        EcoFunction == "P2" | 
                                        EcoFunction == "C1" | 
                                        EcoFunction == "C2" | 
                                        EcoFunction == "Ntot" | 
                                        EcoFunction == "Ptot" | 
                                        EcoFunction == "Ctot", 
                                      "Stock", 
                                      if_else(EcoFunction == "FLUX_P1" | 
                                                EcoFunction == "FLUX_P2" | 
                                                EcoFunction == "FLUX_C1" | 
                                                EcoFunction == "FLUX_C2" | 
                                                EcoFunction == "FLUX_Ptot" | 
                                                EcoFunction == "FLUX_Ctot", 
                                              "Nutrient Flux", 
                                              ifelse(EcoFunction == "CR_Eco1" | 
                                                      EcoFunction == "CR_Eco2" |
                                                      EcoFunction == "CR_MetaEco",
                                                    "C:R",
                                                    "Trophic Productivity")))) %>% 
  dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | 
                                     EcoFunction == "P1" | 
                                     EcoFunction == "C1" | 
                                     EcoFunction == "FLUX_P1" | 
                                     EcoFunction == "FLUX_C1" | 
                                     EcoFunction == "PROD_P1" | 
                                     EcoFunction == "PROD_C1" |
                                     EcoFunction == "CR_Eco1", 
                                   "Ecosystem 1", 
                                   if_else(EcoFunction == "N2" | 
                                             EcoFunction == "P2" | 
                                             EcoFunction == "C2" | 
                                             EcoFunction == "FLUX_P2" | 
                                             EcoFunction == "FLUX_C2" | 
                                             EcoFunction == "PROD_P2" | 
                                             EcoFunction == "PROD_C2" |
                                             EcoFunction == "CR_Eco2", 
                                           "Ecosystem 2", 
                                           "Meta-ecosystem"))) %>% 
  dplyr::mutate(., Fertility = factor(Fertility), 
                Function = factor(Function), 
                Scale = factor(Scale), 
                EcoFunction = factor(EcoFunction))

head(EnvFert_CRbelow1_long)

# A tibble: 6 × 9
  biosense Fertility  TIME    I1    I2 EcoFunction  Value Function
  <fct>    <fct>     <dbl> <dbl> <dbl> <fct>        <dbl> <fct>   
1 yes      Equal         2    10    10 N1           0.773 Stock   
2 yes      Equal         2    10    10 P1          16.7   Stock   
3 yes      Equal         2    10    10 C1           0.429 Stock   
4 yes      Equal         2    10    10 N2           0.858 Stock   
5 yes      Equal         2    10    10 P2          12.4   Stock   
6 yes      Equal         2    10    10 C2           3.20  Stock   
# … with 1 more variable: Scale <fct>

And now we use the newly created dataframe EnvFert_CRbelow1_long to produce a figure analogous to the Figure S1 in the Supplementary Materials.

Show code

StockSummCRbelow1 <- EnvFert_CRbelow1_long %>% 
  filter(., Function == "Stock") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, 
                                            Nutrients = "N1", 
                                            Nutrients = "N2", 
                                            Nutrients = "Ntot", 
                                            "Primary\n Producers" = "P1", 
                                            "Primary\n Producers" = "P2", 
                                            "Primary\n Producers" = "Ptot", 
                                            Consumers = "C1", 
                                            Consumers = "C2", 
                                            Consumers = "Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, 
                                             "Nutrients", 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 10, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) +
  scale_color_highcontrast(reverse = F, 
                           name = "Consumer movement", 
                           labels = c("Equal", 
                                      "Along-gradient", 
                                      "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Stock and Biomass") +
  ylab(" ") +
  xlab(" ") +
  coord_cartesian(ylim = c(0, 75)) 

FluxSummCRbelow1 <- EnvFert_CRbelow1_long %>% 
  filter(., Function == "Nutrient Flux") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, 
                                            "Primary\n Producers" = "FLUX_P1", 
                                            "Primary\n Producers" = "FLUX_P2", 
                                            "Primary\n Producers" = "FLUX_Ptot", 
                                            Consumers = "FLUX_C1", 
                                            Consumers = "FLUX_C2", 
                                            Consumers = "FLUX_Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Equal", 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) + 
  scale_color_highcontrast(reverse = F, 
                           name = "Consumer movement", 
                           labels = c("Equal", 
                                      "Along-gradient", 
                                      "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Flux") +
  ylab(" ") +
  xlab(" ") +
  coord_cartesian(ylim = c(0, 350))

ProdSummCRbelow1 <- EnvFert_CRbelow1_long %>% 
  filter(., Function == "Trophic Productivity") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, 
                                            "Primary\n Producers" = "PROD_P1", 
                                            "Primary\n Producers" = "PROD_P2", 
                                            "Primary\n Producers" = "PROD_Ptot", 
                                            Consumers = "PROD_C1", 
                                            Consumers = "PROD_C2", 
                                            Consumers = "PROD_Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Equal", 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) + 
  scale_color_highcontrast(reverse = F, 
                           name = "Consumer movement", 
                           labels = c("Equal", 
                                      "Along-gradient", 
                                      "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Trophic Compartment Productivity") +
  ylab(" ") +
  xlab("Trophic Compartment") +
  coord_cartesian(ylim = c(0, 500)) 

FuncAllCRbelow1 <- StockSummCRbelow1 + FluxSummCRbelow1 + ProdSummCRbelow1 +
  plot_layout(ncol = 1, nrow = 3, guides = "collect") + 
  plot_annotation( tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
  consmov_theme

# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering parameter sets that produce Consumer:Resource ratios below 1.",

FuncAllCRbelow1

Figure 20: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameters sets for which C:R < 1 (unstransformed data). Notably, we observe a clear quantative difference in the nutrient stock at the meta-ecosystem scale, compared with Figure S1. However, this does not appear to cause a qualitative change in the patterns of local and meta-ecosystem functioning observe, as noted in Figure S2. All specifications as in Figure S1.

Show code

ggsave(FuncAllCRbelow1, filename = "../Results/FunctionSummary_CRbelow1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)

While the quantitative changes are more noticeable here, again in the distance between the upper and lower hinges of the boxes, there does not seem to be a marked change in the qualitative pattern observed in our results for any of the ecosystem functions of interest.

Let’s proceed to work on the \(LRR_X\) values.

Log-transformed Response Ratio

As we just did for the untransformed data, let’s use the previously created object EnvFertRR to select the parameter sets producing \(C{:}R < 1\) and store them in a new object, EnvFertRR_CRbelow1.

Show code

EnvFertRR_CRbelow1 <- dplyr::filter(EnvFertRR, CR_Eco1 <= 1 & CR_Eco2 <= 1 & CR_MetaEco <= 1)

We now lengthen this newly created object.

Show code

EnvFertRR_CRbelow1_long <- select(EnvFertRR_CRbelow1, !(N1:e2)) %>% 
  pivot_longer(., 
               cols = N1rr:CR_MetaEco, 
               names_to = "Compartment",
               values_to = "Ratio") %>% 
  dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | 
                                        Compartment == "FLUX_P2rr" | 
                                        Compartment == "FLUX_C1rr" | 
                                        Compartment == "FLUX_C2rr" | 
                                        Compartment == "FLUX_Prr" | 
                                        Compartment == "FLUX_Crr", 
                                      "Nutrient Flux", 
                                      if_else(Compartment == "PROD_P1rr" | 
                                                Compartment == "PROD_C1rr" | 
                                                Compartment == "PROD_P2rr" | 
                                                Compartment == "PROD_C2rr" | 
                                                Compartment == "PROD_Prr" | 
                                                Compartment == "PROD_Crr", 
                                              "Trophic Productivity", 
                                              ifelse(Compartment == "CR_Eco1" | 
                                                      Compartment == "CR_Eco2" |
                                                      Compartment == "CR_MetaEco",
                                                    "C:R",
                                                    "Stock")))) %>% 
  dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | 
                                     Compartment == "P1rr" | 
                                     Compartment == "C1rr" | 
                                     Compartment == "FLUX_P1rr" | 
                                     Compartment == "FLUX_C1rr" | 
                                     Compartment == "PROD_P1rr" | 
                                     Compartment == "PROD_C1rr" |
                                     Compartment == "CR_Eco1", 
                                   "Ecosystem 1", 
                                   if_else(Compartment == "N2rr" | 
                                             Compartment == "P2rr" | 
                                             Compartment == "C2rr" | 
                                             Compartment == "FLUX_P2rr" | 
                                             Compartment == "FLUX_C2rr" | 
                                             Compartment == "PROD_P2rr" | 
                                             Compartment == "PROD_C2rr" |
                                             Compartment == "CR_Eco2", 
                                           "Ecosystem 2", 
                                           "Meta-ecosystem"))) %>% 
  dplyr::mutate(., Fertility = factor(Fertility), 
                Function = factor(Function), 
                Scale = factor(Scale), 
                Compartment = factor(Compartment)) %>%
  dplyr::mutate(., biosense = fct_drop(biosense), 
                stable = fct_drop(stable))

And finally we produce Figure S2, which is shown in the Supplementary Materials. Accordingly, we show the code here but not the figure itself.

Show code

# we make sure the y-axis of this figure has the same range of the "all data"
# graph above by extracting the limits of the y-axis from that graph and
# using them in this one
LRRstock_ymin <- layer_scales(EFstockLRR)$y$get_limits()[1]
LRRstock_ymax <- layer_scales(EFstockLRR)$y$get_limits()[2]

EFstockLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>% 
  filter(., Function == "Stock") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, 
                                            Nutrients = "N1rr", 
                                            Nutrients = "N2rr", 
                                            Nutrients = "Nrr", 
                                            "Primary\n Producers" = "P1rr", 
                                            "Primary\n Producers" = "P2rr", 
                                            "Primary\n Producers" = "Prr", 
                                            Consumers = "C1rr", 
                                            Consumers = "C2rr", 
                                            Consumers = "Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, 
                                             "Nutrients", 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), 
            fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, 
                     name = "Consumer movement", 
                     labels = c("Along-gradient", "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Stock and Biomass") +
  ylab("LRR") +
  xlab(" ") +
  coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))

# again, we grab the y-axis limits from the relevant all-data graph
LRRflux_ymin <- layer_scales(EFfluxLRR)$y$get_limits()[1]
LRRflux_ymax <- layer_scales(EFfluxLRR)$y$get_limits()[2]

EFfluxLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>% 
  filter(., Function == "Nutrient Flux") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, 
                                            "Primary\n Producers" = "FLUX_P1rr", 
                                            "Primary\n Producers" = "FLUX_P2rr", 
                                            "Primary\n Producers" = "FLUX_Prr", 
                                            Consumers = "FLUX_C1rr", 
                                            Consumers = "FLUX_C2rr", 
                                            Consumers = "FLUX_Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), 
            fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, 
                     name = "Consumer movement", 
                     labels = c("Along-gradient", 
                                "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Flux") +
  ylab("LRR") +
  xlab(" ") +
  coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))

# again, we grab the y-axis limits from the relevant all-data graph
LRRprod_ymin <- layer_scales(EFprodLRR)$y$get_limits()[1]
LRRprod_ymax <- layer_scales(EFprodLRR)$y$get_limits()[2]

EFprodLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>% 
  filter(., Function == "Trophic Productivity") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, 
                                            "Primary\n Producers" = "PROD_P1rr", 
                                            "Primary\n Producers" = "PROD_P2rr", 
                                            "Primary\n Producers" = "PROD_Prr", 
                                            Consumers = "PROD_C1rr", 
                                            Consumers = "PROD_C2rr", 
                                            Consumers = "PROD_Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), 
            fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) +
  scale_color_manual(values = ptHCsubset, 
                    name = "Consumer movement", 
                    labels = c("Along-gradient", 
                               "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Trophic Compartment Productivity") +
  ylab("LRR") +
  xlab("Trophic compartment") + 
  coord_cartesian(ylim = c(LRRprod_ymin, LRRprod_ymax))

EFallLRR_CRbelow1 <- EFstockLRR_CRbelow1 + EFfluxLRR_CRbelow1 + EFprodLRR_CRbelow1 + 
  plot_layout(ncol = 1, nrow = 3, guides = "collect") + 
  plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
  consmov_theme

# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering the LRR of parameter sets that produce Consumer:Resource ratios below 1.", 

# EFallLRR_CRbelow1

ggsave(EFallLRR_CRbelow1, filename = "../Results/FunctionSummary_LRR_CRbelow1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)

And here, we produce Figure S4 from the Supplementary Materials, showing median parameter values for the subset of results where \(C{:}R < 1\). Again, we show the code but not the figure itself.

Show code

# table of median parameter values
# EnvFertRR_CRbelow1 %>% select(., Fertility, u1:e2) %>% ungroup() %>% group_by(., Fertility) %>% pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>% mutate(., Parameter = as_factor(Parameter)) %>% ungroup() %>% group_by(., Fertility, Parameter) %>% summarise(., across(.cols = value, .fns = list(mean = mean, sd = sd)), .groups = "keep") %>% pivot_wider(., names_from = "Fertility", values_from = c("value_mean", "value_sd"))%>% gt(rowname_col = "Parameter", groupname_col = "Fertility") %>% fmt_number(., columns = 2:5, decimals = 2) %>% tab_spanner(label = md("Along"), columns = c("value_mean_Along", "value_sd_Along")) %>% tab_spanner(label = md("Against"), columns = c("value_mean_Against", "value_sd_Against")) %>% cols_label(value_mean_Against = "mean", value_sd_Against = "sd", value_mean_Along = "mean", value_sd_Along = "sd") 

# %>% gtsave(., filename = "../Results/MeanParam_LRR_CRbelow1.rtf")

EnvFertRR_CRbelow1 %>% 
  select(., Fertility, u1:e2) %>% 
  ungroup() %>% 
  group_by(., Fertility) %>% 
  pivot_longer(., cols = -Fertility, 
               names_to = "Parameter", 
               values_to = "value") %>% 
  mutate(., Parameter = as_factor(Parameter)) %>% 
  ggplot(., aes(x = Parameter, y = value)) +
  geom_boxjitter(aes(fill = Parameter, 
                     col = Parameter), 
                 jitter.fill = NULL, alpha = 0.1) + 
  scale_color_manual(values=met.brewer("Signac", 13)) + 
  scale_fill_manual(values=met.brewer("Signac", 13)) + 
  theme_pubr(legend = "none") + 
  ylab("Mean") + 
  facet_grid(.~Fertility, 
             labeller = as_labeller(c(`Against` = "Against-gradient Consumer movement", 
                                      `Along` = "Along-gradient Consumer movement")))  

Show code

# ggsave(filename = "../Results/MeanParam_LRR_CRbelow1.pdf", dpi = 300, device = "pdf")

Consumer movement in top-heavy systems

Here we reproduce Figures 2, 3, and 4 from the main text and Figure S1 from the Supplementary Materials using only parameter sets that produce \(C{:}R \in (1,10)\). We exclude \(C{:}R > 10\) values, as these are likely extremely rare in real-world ecosystems. We begin, as usual, with the untransformed dataframe.

Untransformed parameter sets

First, let’s exclude from EnvFert_respos all the parameter sets producing \(C{:}R > 10\) and those producing \(C{:}R < 1\).

Show code

EnvFert_CRabove1 <- dplyr::mutate(EnvFert_respos,
                         CR_Eco1 = C1/P1, # C:R for Ecosystem 1
                         CR_Eco2 = C2/P2, # C:R for Ecosystem 2
                         CR_MetaEco = Ctot/Ptot # C:R for meta-ecosystem
                         ) %>%
  dplyr::filter(., CR_Eco1 <= 10 & CR_Eco2 <= 10 & CR_MetaEco <= 10) %>% 
  dplyr::filter(., CR_Eco1 > 1 & CR_Eco2 > 1 & CR_MetaEco > 1)

Now we lengthen the dataset EnvFert_CRabove1 created in the code chunk above for plotting.

Show code

EnvFert_CRabove1_long <- select(EnvFert_CRabove1, !u1:e2) %>% 
  pivot_longer(., 
               cols = N1:CR_MetaEco, 
               names_to = "EcoFunction", 
               values_to = "Value") %>% 
  dplyr::mutate(., Function = if_else(EcoFunction == "N1" | 
                                        EcoFunction == "N2" | 
                                        EcoFunction == "P1" | 
                                        EcoFunction == "P2" | 
                                        EcoFunction == "C1" | 
                                        EcoFunction == "C2" | 
                                        EcoFunction == "Ntot" | 
                                        EcoFunction == "Ptot" | 
                                        EcoFunction == "Ctot", 
                                      "Stock", 
                                      if_else(EcoFunction == "FLUX_P1" | 
                                                EcoFunction == "FLUX_P2" | 
                                                EcoFunction == "FLUX_C1" | 
                                                EcoFunction == "FLUX_C2" | 
                                                EcoFunction == "FLUX_Ptot" | 
                                                EcoFunction == "FLUX_Ctot", 
                                              "Nutrient Flux", 
                                              ifelse(EcoFunction == "CR_Eco1" | 
                                                      EcoFunction == "CR_Eco2" |
                                                      EcoFunction == "CR_MetaEco",
                                                    "C:R",
                                                    "Trophic Productivity")))) %>% 
  dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | 
                                     EcoFunction == "P1" | 
                                     EcoFunction == "C1" | 
                                     EcoFunction == "FLUX_P1" | 
                                     EcoFunction == "FLUX_C1" | 
                                     EcoFunction == "PROD_P1" | 
                                     EcoFunction == "PROD_C1" |
                                     EcoFunction == "CR_Eco1", 
                                   "Ecosystem 1", 
                                   if_else(EcoFunction == "N2" | 
                                             EcoFunction == "P2" | 
                                             EcoFunction == "C2" | 
                                             EcoFunction == "FLUX_P2" | 
                                             EcoFunction == "FLUX_C2" | 
                                             EcoFunction == "PROD_P2" | 
                                             EcoFunction == "PROD_C2" |
                                             EcoFunction == "CR_Eco2", 
                                           "Ecosystem 2", 
                                           "Meta-ecosystem"))) %>% 
  dplyr::mutate(., Fertility = factor(Fertility), 
                Function = factor(Function), 
                Scale = factor(Scale), 
                EcoFunction = factor(EcoFunction))

head(EnvFert_CRabove1_long)

# A tibble: 6 × 9
  biosense Fertility  TIME    I1    I2 EcoFunction  Value Function
  <fct>    <fct>     <dbl> <dbl> <dbl> <fct>        <dbl> <fct>   
1 yes      Equal        13    10    10 N1           7.29  Stock   
2 yes      Equal        13    10    10 P1           0.848 Stock   
3 yes      Equal        13    10    10 C1           3.26  Stock   
4 yes      Equal        13    10    10 N2          15.5   Stock   
5 yes      Equal        13    10    10 P2           0.894 Stock   
6 yes      Equal        13    10    10 C2           7.23  Stock   
# … with 1 more variable: Scale <fct>

And, finally, we reproduce Figure S1 from the Supplementary Materials.

Show code

StockSummCRabove1 <- EnvFert_CRabove1_long %>% 
  filter(., Function == "Stock") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, 
                                            Nutrients = "N1", 
                                            Nutrients = "N2", 
                                            Nutrients = "Ntot", 
                                            "Primary\n Producers" = "P1", 
                                            "Primary\n Producers" = "P2", 
                                            "Primary\n Producers" = "Ptot", 
                                            Consumers = "C1", 
                                            Consumers = "C2", 
                                            Consumers = "Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, 
                                             "Nutrients", 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along", "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 10, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) +
  scale_color_highcontrast(reverse = F, 
                           name = "Consumer movement", 
                           labels = c("Equal", 
                                      "Along-gradient", 
                                      "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Stock and Biomass") +
  ylab(" ") +
  xlab(" ") +
  coord_cartesian(ylim = c(0, 75)) 

FluxSummCRabove1 <- EnvFert_CRabove1_long %>% 
  filter(., Function == "Nutrient Flux") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, 
                                            "Primary\n Producers" = "FLUX_P1", 
                                            "Primary\n Producers" = "FLUX_P2", 
                                            "Primary\n Producers" = "FLUX_Ptot", 
                                            Consumers = "FLUX_C1", 
                                            Consumers = "FLUX_C2", 
                                            Consumers = "FLUX_Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Equal", 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) + 
  scale_color_highcontrast(reverse = F, 
                           name = "Consumer movement", 
                           labels = c("Equal", 
                                      "Along-gradient", 
                                      "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Flux") +
  ylab(" ") +
  xlab(" ") +
  coord_cartesian(ylim = c(0, 350))

ProdSummCRabove1 <- EnvFert_CRabove1_long %>% 
  filter(., Function == "Trophic Productivity") %>% 
  dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, 
                                            "Primary\n Producers" = "PROD_P1", 
                                            "Primary\n Producers" = "PROD_P2", 
                                            "Primary\n Producers" = "PROD_Ptot", 
                                            Consumers = "PROD_C1", 
                                            Consumers = "PROD_C2", 
                                            Consumers = "PROD_Ctot")) %>%
  dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Equal", 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = EcoFunction, y = Value)) +
  geom_hline(yintercept = 50, linetype = "dashed", lwd = 0.25) + 
  geom_boxplot(aes(color = Fertility)) + 
  scale_color_highcontrast(reverse = F, 
                           name = "Consumer movement", 
                           labels = c("Equal", 
                                      "Along-gradient", 
                                      "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Trophic Compartment Productivity") +
  ylab(" ") +
  xlab("Trophic compartment") +
  coord_cartesian(ylim = c(0, 500)) 

FuncAllCRabove1 <- StockSummCRabove1 + FluxSummCRabove1 + ProdSummCRabove1 + 
  plot_layout(ncol = 1, nrow = 3, guides = "collect") + 
  plot_annotation( tag_levels = "a", tag_prefix = "(", tag_suffix = ")") & 
  consmov_theme

# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering parameter sets that produce Consumer:Resource ratios above 1.",

FuncAllCRabove1

Figure 21: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameter sets for which C:R > 1 and < 10 (unstransformed data). Here, biomass and nutrient stock results appear to more clearly point towards a trophic cascade at play at the meta-ecosystem scale, compared to Figure S1. Patterns of nutrient flux and trophic compartment productivity appear qualitatively unchanged. All specifications as in Figure S1.

Show code

ggsave(FuncAllCRabove1, filename = "../Results/FunctionSummary_CRabove1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)

Compared with Figure S1, we see a clear change in the pattern of biomass accumulation—consistent with the inverted, top-heavy biomass pyramid expected when \(C{:}R > 1\) (McCauley et al. 2018).

Log-transformed Response Ratio

Here, we produce Figure S3 from the Supplementary Materials. We use only parameter sets that produce \(C{:}R \in (1,10)\). First, we exclude all \(C{:}R\) values below 1.

Show code

EnvFertRR_CRabove1 <- dplyr::filter(EnvFertRR, CR_Eco1 <= 10 & CR_Eco2 <= 10 & CR_MetaEco <= 10) %>% dplyr::filter(., CR_Eco1 > 1 & CR_Eco2 > 1 & CR_MetaEco > 1)

Then, we lengthen the dataset for plotting.

Show code

EnvFertRR_CRabove1_long <- select(EnvFertRR_CRabove1, !(N1:e2)) %>% 
  pivot_longer(., 
               cols = N1rr:CR_MetaEco, 
               names_to = "Compartment",
               values_to = "Ratio") %>% 
  dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | 
                                        Compartment == "FLUX_P2rr" | 
                                        Compartment == "FLUX_C1rr" | 
                                        Compartment == "FLUX_C2rr" | 
                                        Compartment == "FLUX_Prr" | 
                                        Compartment == "FLUX_Crr", 
                                      "Nutrient Flux", 
                                      if_else(Compartment == "PROD_P1rr" | 
                                                Compartment == "PROD_C1rr" | 
                                                Compartment == "PROD_P2rr" | 
                                                Compartment == "PROD_C2rr" | 
                                                Compartment == "PROD_Prr" | 
                                                Compartment == "PROD_Crr", 
                                              "Trophic Productivity", 
                                              ifelse(Compartment == "CR_Eco1" | 
                                                      Compartment == "CR_Eco2" |
                                                      Compartment == "CR_MetaEco",
                                                    "C:R",
                                                    "Stock")))) %>% 
  dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | 
                                     Compartment == "P1rr" | 
                                     Compartment == "C1rr" | 
                                     Compartment == "FLUX_P1rr" | 
                                     Compartment == "FLUX_C1rr" | 
                                     Compartment == "PROD_P1rr" | 
                                     Compartment == "PROD_C1rr" |
                                     Compartment == "CR_Eco1", 
                                   "Ecosystem 1", 
                                   if_else(Compartment == "N2rr" | 
                                             Compartment == "P2rr" | 
                                             Compartment == "C2rr" | 
                                             Compartment == "FLUX_P2rr" | 
                                             Compartment == "FLUX_C2rr" | 
                                             Compartment == "PROD_P2rr" | 
                                             Compartment == "PROD_C2rr" |
                                             Compartment == "CR_Eco2", 
                                           "Ecosystem 2", 
                                           "Meta-ecosystem"))) %>% 
  dplyr::mutate(., Fertility = factor(Fertility), 
                Function = factor(Function), 
                Scale = factor(Scale), 
                Compartment = factor(Compartment)) %>%
  dplyr::mutate(., biosense = fct_drop(biosense), 
                stable = fct_drop(stable))

And now we produce Figure S3. Note that, as above, we do not show the figure itself here but only the code needed to produce it.

Show code

EFstockLRR_CRabove1 <- EnvFertRR_CRabove1_long %>% 
  filter(., Function == "Stock") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, 
                                            Nutrients = "N1rr", 
                                            Nutrients = "N2rr", 
                                            Nutrients = "Nrr", 
                                            "Primary\n Producers" = "P1rr", 
                                            "Primary\n Producers" = "P2rr", 
                                            "Primary\n Producers" = "Prr", 
                                            Consumers = "C1rr", 
                                            Consumers = "C2rr", 
                                            Consumers = "Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, 
                                             "Nutrients", 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), 
            fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, 
                     name = "Consumer movement", 
                     labels = c("Along-gradient", "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Stock and Biomass") +
  ylab("LRR") +
  xlab(" ") +
  coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))

EFfluxLRR_CRabove1 <- EnvFertRR_CRabove1_long %>% 
  filter(., Function == "Nutrient Flux") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, 
                                            "Primary\n Producers" = "FLUX_P1rr", 
                                            "Primary\n Producers" = "FLUX_P2rr", 
                                            "Primary\n Producers" = "FLUX_Prr", 
                                            Consumers = "FLUX_C1rr", 
                                            Consumers = "FLUX_C2rr", 
                                            Consumers = "FLUX_Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), 
            fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, 
                     name = "Consumer movement", 
                     labels = c("Along-gradient", 
                                "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Nutrient Flux") +
  ylab("LRR") +
  xlab(" ") +
  coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))

EFprodLRR_CRabove1 <- EnvFertRR_CRabove1_long %>% 
  filter(., Function == "Trophic Productivity") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, 
                                            "Primary\n Producers" = "PROD_P1rr", 
                                            "Primary\n Producers" = "PROD_P2rr", 
                                            "Primary\n Producers" = "PROD_Prr", 
                                            Consumers = "PROD_C1rr", 
                                            Consumers = "PROD_C2rr", 
                                            Consumers = "PROD_Crr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, 
                                             "Primary\n Producers", 
                                             "Consumers")) %>%
  dplyr::mutate(., Fertility = fct_relevel(Fertility, 
                                           "Along", 
                                           "Against")) %>%
  ggplot(., aes(x = Compartment, y = log10(Ratio))) +
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0), 
            fill = "gray95") +
  geom_hline(yintercept = 0, linetype = "dashed", colour = "black") + 
  geom_boxplot(aes(color = Fertility, fill = Fertility)) + 
  scale_color_manual(values = ptHCsubset, 
                    name = "Consumer movement", 
                    labels = c("Along-gradient", 
                               "Against-gradient")) +
  scale_fill_manual(values = ptHCsubset_fill, name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
  facet_grid(.~Scale, scales = "free") +
  labs(title = "Trophic Compartment Productivity") +
  ylab("LRR") +
  xlab("Trophic compartment") + 
  coord_cartesian(ylim = c(LRRprod_ymin, LRRprod_ymax))

EFallLRR_CRabove1 <- EFstockLRR_CRabove1 + EFfluxLRR_CRabove1 + EFprodLRR_CRabove1 + 
  plot_layout(ncol = 1, nrow = 3, guides = "collect") + 
  plot_annotation(tag_levels = "a", tag_prefix = "(", tag_suffix = ")") &
  consmov_theme

# title = "Consumer movement influence on local and meta-ecosystem functions", subtitle = "When considering the LRR of parameter sets that produce Consumer:Resource ratios above 1 but below 10.",

# EFallLRR_CRabove1

ggsave(EFallLRR_CRabove1, filename = "../Results/FunctionSummary_LRR_CRabove1.pdf", device = "pdf", dpi = 300, width = 10, height = 10)

As above, we now plot the median parameter values for these parameter sets, which are shown in Figure S5 in the Supplementary Materials. Accordingly, we only show the code to produce the figure and not the figure itself.

Show code

# table of median parameter sets

# EnvFertRR_CRabove1 %>% select(., Fertility, u1:e2) %>% ungroup() %>% group_by(., Fertility) %>% pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>% mutate(., Parameter = as_factor(Parameter)) %>% ungroup() %>% group_by(., Fertility, Parameter) %>% summarise(., across(.cols = value, .fns = list(mean = mean, sd = sd)), .groups = "keep") %>% pivot_wider(., names_from = "Fertility", values_from = c("value_mean", "value_sd"))%>% gt(rowname_col = "Parameter", groupname_col = "Fertility") %>% fmt_number(., columns = 2:5, decimals = 2) %>% tab_spanner(label = md("Along"), columns = c("value_mean_Along", "value_sd_Along")) %>% tab_spanner(label = md("Against"), columns = c("value_mean_Against", "value_sd_Against")) %>% cols_label(value_mean_Against = "mean", value_sd_Against = "sd", value_mean_Along = "mean", value_sd_Along = "sd") 
EnvFertRR_CRabove1 %>% 
  select(., Fertility, u1:e2) %>% 
  ungroup() %>%
  group_by(., Fertility) %>% 
  pivot_longer(., cols = -Fertility, 
               names_to = "Parameter", 
               values_to = "value") %>% 
  mutate(., Parameter = as_factor(Parameter)) %>% 
  ggplot(., aes(x = Parameter, y = value)) + 
  geom_boxjitter(aes(fill = Parameter, col = Parameter), 
                 jitter.fill = NULL, alpha = 0.1) + 
  scale_color_manual(values=met.brewer("Signac", 13)) + 
  scale_fill_manual(values=met.brewer("Signac", 13)) + 
  theme_pubr(legend = "none") +
  ylab("Mean") + 
  facet_grid(.~Fertility, 
             labeller = as_labeller(c(`Against` = "Against-gradient Consumer movement", 
                                      `Along` = "Along-gradient Consumer movement")))  

Show code

# ggsave(filename = "../Results/MeanParam_LRR_CRabove1.pdf", dpi = 300, device = "pdf")

Summary tables

Here, we produce summary tables for the results of the model’s iterations, using the \(LRR_X\) and \(C{:}R\) values we calculated to produce the summary plots above. For \(LRR_X\), we report median values for each ecosystem function of interest in both local ecosystems and in the meta-ecosystem.

For ease of interpretation, median \(LRR_X\) values \(< 0\) will be shown in shades of magenta based on how close they are to -1 and median \(LRR_X > 0\) will be shown in shades of green based on their closeness to 1.

First, we produce dataframes to contain the median values, one for each of stock, nutrient flux, and primary and secondary productivity.

Show code

# first, compute median of log10 ratio values

EnvFertLRR_tables <- dplyr::filter(EnvFertRR_long, Ratio>=0) %>% dplyr::group_by(., Fertility, Function, Scale, Compartment) %>% dplyr::mutate(., LRR = log10(Ratio)) %>% dplyr::summarise(., min = min(LRR), median = median(LRR), max = max(LRR), .groups = "keep")

Now, we produce the tables, separating between along and against-gradient movement, as above.

Consumer movement changes local and meta-ecosystem functions

Biomass and Nutrient Stock

Show code

EFStockTab <- EnvFertLRR_tables %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(., Function == "Stock") %>% 
  pivot_wider(., names_from = Fertility, values_from = c(min, median, max)) %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "C1rr", "Consumers" = "C2rr", "Consumers" = "Crr", "Producers" = "P1rr", "Producers" = "P2rr", "Producers" = "Prr", "Nutrients" = "N1rr", "Nutrients" = "N2rr", "Nutrients" = "Nrr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients", "Producers", "Consumers")) %>%
  gt(groupname_col = "Scale") %>% 
  cols_hide("Function") %>% 
  tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along", "median_Along", "max_Along")) %>%
  tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against", "median_Against", "max_Against")) %>%
  cols_label(Compartment = md("**Compartment**"), "min_Along" = md("*Min.*"), "median_Along" = md("*Median*"), "max_Along" = md("*Max.*"), "min_Against" = md("*Min.*"), "median_Against" = md("*Median*"), "max_Against" = md("*Max.*")) %>%
  fmt_number(., columns = c(4:9), decimals = 2) %>% 
  data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>% tab_header(title = md("**Change in nutrient stock and organic biomass accumulation** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."), subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>% tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

# gtsave(EFStockTab, filename = "EFbiostock.tex", path = "../Results/")

EFStockTab

Change in nutrient stock and organic biomass accumulation in the meta-ecosystem when consumers move against or along an environmental fertility gradient.
Values are expressed as median LRR for 10000 iterations of the model.
Compartment	Against-gradient			Along-gradient
	Min.	Median	Max.	Min.	Median	Max.
Ecosystem 1
Consumers	−3.16	−0.72	−0.70	0.26	0.26	2.45
Nutrients	−0.70	−0.16	0.00	0.00	0.12	0.26
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	−0.80	0.17	0.92	−1.35	−0.27	0.31
Nutrients	−0.68	0.11	0.26	−0.70	−0.14	0.25
Producers	0.00	0.04	1.46	−2.72	−0.12	0.00
Meta-ecosystem
Consumers	−1.16	0.08	0.59	−1.08	−0.09	0.57
Nutrients	−0.67	0.03	0.25	−0.68	−0.03	0.25
Producers	0.00	0.01	0.90	−0.80	−0.02	0.00

Nutrient Flux

Show code

EFFluxTab <- EnvFertLRR_tables %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(., Function == "Nutrient Flux") %>% 
  pivot_wider(., names_from = Fertility, values_from = c(min, median, max)) %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "FLUX_C1rr", "Consumers" = "FLUX_C2rr", "Consumers" = "FLUX_Crr", "Producers" = "FLUX_P1rr", "Producers" = "FLUX_P2rr", "Producers" = "FLUX_Prr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
  gt(groupname_col = "Scale") %>% 
  cols_hide("Function") %>% 
  tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along", "median_Along", "max_Along")) %>%
  tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against", "median_Against", "max_Against")) %>%
  cols_label(Compartment = md("**Compartment**"), "min_Along" = md("*Min.*"), "median_Along" = md("*Median*"), "max_Along" = md("*Max.*"), "min_Against" = md("*Min.*"), "median_Against" = md("*Median*"), "max_Against" = md("*Max.*"))%>%
  fmt_number(., columns = c(4:9), decimals = 2) %>% 
  data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>% tab_header(title = md("**Change in nutrient flux** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."), subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>% tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

# gtsave(EFFluxTab, filename = "EFnutflux.tex", path = "../Results/")

EFFluxTab

Change in nutrient flux in the meta-ecosystem when consumers move against or along an environmental fertility gradient.
Values are expressed as median LRR for 10000 iterations of the model.
Compartment	Against-gradient			Along-gradient
	Min.	Median	Max.	Min.	Median	Max.
Ecosystem 1
Consumers	−3.16	−0.72	−0.70	0.26	0.26	2.45
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	−0.80	0.17	0.92	−1.35	−0.27	0.31
Producers	0.00	0.04	1.46	−2.72	−0.12	0.00
Meta-ecosystem
Consumers	−0.86	0.07	0.49	−1.06	−0.09	0.57
Producers	0.00	0.01	0.90	−0.86	−0.02	0.00

Trophic Productivity

Show code

EFProdTab <- EnvFertLRR_tables %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(., Function == "Trophic Productivity") %>% 
  pivot_wider(., names_from = Fertility, values_from = c(min, median, max)) %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "PROD_C1rr", "Consumers" = "PROD_C2rr", "Consumers" = "PROD_Crr", "Producers" = "PROD_P1rr", "Producers" = "PROD_P2rr", "Producers" = "PROD_Prr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
  gt(groupname_col = "Scale") %>% 
  cols_hide("Function") %>% 
  tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along", "median_Along", "max_Along")) %>%
  tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against", "median_Against", "max_Against")) %>%
  cols_label(Compartment = md("**Compartment**"), "min_Along" = md("*Min.*"), "median_Along" = md("*Median*"), "max_Along" = md("*Max.*"), "min_Against" = md("*Min.*"), "median_Against" = md("*Median*"), "max_Against" = md("*Max.*"))%>%
  fmt_number(., columns = c(4:9), decimals = 2) %>% 
  data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>% tab_header(title = md("**Change in trophic compartment productivity** in the meta-ecosystem when consumers move _along-_ or _against-gradient_ of nutrient availability."), subtitle = md("Values are expressed as median _LRR_ for 10 000 iterations of the model.")) %>% tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())

# gtsave(EFProdTab, filename = "EFprod.tex", path = "../Results/")

EFProdTab

Change in trophic compartment productivity in the meta-ecosystem when consumers move along- or against-gradient of nutrient availability.
Values are expressed as median LRR for 10 000 iterations of the model.
Compartment	Against-gradient			Along-gradient
	Min.	Median	Max.	Min.	Median	Max.
Ecosystem 1
Consumers	−3.16	−0.72	−0.70	0.26	0.26	2.45
Producers	−0.70	−0.16	0.00	0.00	0.12	0.26
Ecosystem 2
Consumers	0.01	0.22	1.40	−2.63	−0.47	−0.01
Producers	0.00	0.18	1.40	−2.63	−0.34	0.00
Meta-ecosystem
Consumers	−0.88	0.03	0.47	−1.02	−0.04	0.70
Producers	−0.47	0.03	0.88	−0.81	−0.04	0.22

Consumer movement changes biomass distribution at multiple spatial scales

Movement from ecosystem 1

Show code

# first, compute median of log10 ratio values

ecoCR_g_table <- bind_rows("above 1" = ecoCRg_lmOutput_above1, 
                           "below 1" = ecoCRg_lmOutput_below1, 
                           .id = "CRslice") %>% 
  ungroup() %>% 
  select(., CRslice:adj.r.squared, g_slope:g_se) %>% 
  pivot_wider(., names_from = "CRslice", 
              values_from = c(r.squared, adj.r.squared, g_slope, g_se)) %>% 
  gt(., groupname_col = "Scale") %>% 
  tab_spanner(label = md("**C:R (1, 10]**"), 
              columns = c("r.squared_above 1", 
                          "adj.r.squared_above 1", 
                          "g_slope_above 1",
                          "g_se_above 1")) %>% 
  tab_spanner(label = md("**C:R [0, 1]**"), 
              columns = c("r.squared_below 1", 
                          "adj.r.squared_below 1", 
                          "g_slope_below 1",
                          "g_se_below 1")) %>% 
  fmt_number(., columns = 3:10, decimals = 3) %>% 
  cols_label(., "r.squared_above 1"= md("R sq."), 
             "adj.r.squared_above 1" = md("Adj. R sq."), 
             "g_slope_above 1" = md("g Coeff."), 
             "g_se_above 1" = md("g SE"),
             "r.squared_below 1"= md("R sq."), 
             "adj.r.squared_below 1" = md("Adj. R sq."), 
             "g_slope_below 1" = md("g Coeff."),
             "g_se_below 1" = md("g SE")) %>% 
  # data_color(., columns = c(3:10), 
  #            colors = col_numeric(palette = "PiYG", domain = c(-1,1))) %>% 
  tab_header(title = md("**Relationships between biomass _C:R_ and movement from ecosystem 1.**"), subtitle = md("The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_  values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>% 
  tab_style(style = cell_text(align = "left", size = 18), 
            locations = cells_title())

# gtsave(ecoCR_g_table, filename = "../Results/CR_g_summary.rtf")

ecoCR_g_table

Relationships between biomass C:R and movement from ecosystem 1.
The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0.
Fertility	C:R (1, 10]				C:R [0, 1]
	R sq.	Adj. R sq.	g Coeff.	g SE	R sq.	Adj. R sq.	g Coeff.	g SE
Ecosystem 1
Against	0.083	0.078	−0.728	0.180	0.119	0.119	−0.017	0.180
Along	0.148	0.147	−0.343	0.022	0.010	0.010	−0.009	0.022
Equal	0.177	0.176	−0.537	0.040	0.039	0.039	−0.016	0.040
Ecosystem 2
Against	0.001	0.001	−0.025	0.016	0.000	0.000	0.001	0.016
Along	0.000	0.000	0.018	0.015	0.006	0.005	0.007	0.015
Equal	0.001	0.001	−0.027	0.015	0.001	0.000	0.003	0.015
Meta-ecosystem
Against	0.034	0.033	−0.078	0.012	0.002	0.001	−0.004	0.012
Along	0.048	0.047	−0.100	0.011	0.002	0.002	−0.004	0.011
Equal	0.049	0.048	−0.106	0.012	0.005	0.005	−0.006	0.012

Movement towards ecosystem 2

Show code

# first, compute median of log10 ratio values

ecoCR_m_table <- bind_rows("above 1" = ecoCRm_lmOutput_above1, 
                           "below 1" = ecoCRm_lmOutput_below1, 
                           .id = "CRslice") %>% 
  ungroup() %>% 
  select(., CRslice:adj.r.squared, m_slope:m_se) %>% 
  pivot_wider(., names_from = "CRslice", 
              values_from = c(r.squared, adj.r.squared, m_slope, m_se)) %>% 
  gt(., groupname_col = "Scale") %>% 
  tab_spanner(label = md("**C:R (1, 10]**"), 
              columns = c("r.squared_above 1", 
                          "adj.r.squared_above 1", 
                          "m_slope_above 1",
                          "m_se_above 1")) %>% 
  tab_spanner(label = md("**C:R [0, 1]**"), 
              columns = c("r.squared_below 1", 
                          "adj.r.squared_below 1", 
                          "m_slope_below 1",
                          "m_se_below 1")) %>% 
  fmt_number(., columns = 3:10, decimals = 3) %>% 
  cols_label(., "r.squared_above 1"= md("R sq."), 
             "adj.r.squared_above 1" = md("Adj. R sq."), 
             "m_slope_above 1" = md("m Coeff."), 
             "m_se_above 1" = md("m SE"),
             "r.squared_below 1"= md("R sq."), 
             "adj.r.squared_below 1" = md("Adj. R sq."), 
             "m_slope_below 1" = md("m Coeff."),
             "m_se_below 1" = md("m SE")) %>% 
  # data_color(., columns = c(3:10), 
  #            colors = col_numeric(palette = "PiYG", domain = c(-,))) %>% 
  tab_header(title = md("**Relationships between biomass _C:R_ and movement into ecosystem 2.**"), subtitle = md("The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_  values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>% 
  tab_style(style = cell_text(align = "left", size = 18), 
            locations = cells_title()) 

# gtsave(ecoCR_m_table, filename = "../Results/CR_m_summary.rtf")

ecoCR_m_table

Relationships between biomass C:R and movement into ecosystem 2.
The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0.
Fertility	C:R (1, 10]				C:R [0, 1]
	R sq.	Adj. R sq.	m Coeff.	m SE	R sq.	Adj. R sq.	m Coeff.	m SE
Ecosystem 1
Against	0.003	−0.002	0.036	0.048	0.000	0.000	0.000	0.001
Along	0.000	0.000	−0.014	0.018	0.001	0.001	−0.002	0.001
Equal	0.002	0.001	−0.031	0.026	0.000	0.000	−0.001	0.001
Ecosystem 2
Against	0.000	0.000	−0.008	0.016	0.000	0.000	0.000	0.002
Along	0.001	0.001	0.026	0.015	0.016	0.016	0.012	0.001
Equal	0.000	0.000	−0.001	0.015	0.006	0.006	0.007	0.001
Meta-ecosystem
Against	0.000	−0.001	0.000	0.013	0.000	0.000	0.000	0.001
Along	0.000	0.000	0.010	0.011	0.000	0.000	0.001	0.001
Equal	0.001	0.000	−0.014	0.013	0.000	0.000	0.001	0.001

Death rate in Dispersers’ pool

Show code

# first, compute median of log10 ratio values

ecoCR_c_table <- bind_rows("above 1" = ecoCRc_lmOutput_above1, 
                           "below 1" = ecoCRc_lmOutput_below1, 
                           .id = "CRslice") %>% 
  ungroup() %>% 
  select(., CRslice:adj.r.squared, c_slope:c_se) %>% 
  pivot_wider(., names_from = "CRslice", 
              values_from = c(r.squared, adj.r.squared, c_slope, c_se)) %>% 
  gt(., groupname_col = "Scale") %>% 
  tab_spanner(label = md("**C:R (1, 10]**"), 
              columns = c("r.squared_above 1", 
                          "adj.r.squared_above 1", 
                          "c_slope_above 1",
                          "c_se_above 1")) %>% 
  tab_spanner(label = md("**C:R [0, 1]**"), 
              columns = c("r.squared_below 1", 
                          "adj.r.squared_below 1", 
                          "c_slope_below 1",
                          "c_se_below 1")) %>% 
  fmt_number(., columns = 3:10, decimals = 3) %>% 
  cols_label(., "r.squared_above 1"= md("R sq."), 
             "adj.r.squared_above 1" = md("Adj. R sq."), 
             "c_slope_above 1" = md("c Coeff."), 
             "c_se_above 1" = md("c SE"),
             "r.squared_below 1"= md("R sq."), 
             "adj.r.squared_below 1" = md("Adj. R sq."), 
             "c_slope_below 1" = md("c Coeff."),
             "c_se_below 1" = md("c SE")) %>% 
  # data_color(., columns = c(3:10), 
  #            colors = col_numeric(palette = "PiYG", domain = c(-,))) %>% 
  tab_header(title = md("**Relationships between biomass _C:R_ and death rate in dispersers' pool.**"), subtitle = md("The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_  values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>% 
  tab_style(style = cell_text(align = "left", size = 18), 
            locations = cells_title()) 

gtsave(ecoCR_c_table, filename = "../Results/CR_c_summary.rtf")

ecoCR_c_table

Relationships between biomass C:R and death rate in dispersers' pool.
The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0.
Fertility	C:R (1, 10]				C:R [0, 1]
	R sq.	Adj. R sq.	c Coeff.	c SE	R sq.	Adj. R sq.	c Coeff.	c SE
Ecosystem 1
Against	0.005	−0.001	0.047	0.050	0.000	0.000	−0.001	0.001
Along	0.000	0.000	−0.012	0.019	0.000	0.000	0.002	0.001
Equal	0.000	−0.001	0.000	0.026	0.000	0.000	0.001	0.001
Ecosystem 2
Against	0.000	0.000	0.010	0.016	0.000	0.000	0.000	0.002
Along	0.000	0.000	−0.004	0.014	0.002	0.002	−0.004	0.001
Equal	0.001	0.001	−0.030	0.015	0.002	0.002	−0.005	0.002
Meta-ecosystem
Against	0.003	0.003	0.025	0.012	0.000	0.000	0.000	0.001
Along	0.000	−0.001	−0.003	0.011	0.001	0.001	−0.003	0.001
Equal	0.000	−0.001	0.004	0.012	0.001	0.001	−0.003	0.001

Summary comparison across C:R values

Here, we produce summary tables that show how the effects of consumer movement on local and meta-ecosystem functions vary as the biomass distribution changes in the system as a result of the consumer movement itself.

We begin by collating together the required dataframes.

Show code

CRcomparison <- bind_rows("all" = EnvFertRR, 
                          "above_1" = EnvFertRR_CRabove1, 
                          "below_1" = EnvFertRR_CRbelow1, 
                          .id = "CR") %>% 
  select(., !(N1:Ctot)) %>% 
  filter(., biosense == "yes" | stable == "stable") %>%
  mutate(., CR = as_factor(CR)) %>%
  droplevels()

Then, we lengthen the dataset as this is the preferred format for producing tables with package gt.

Show code

CRcomparison_long <- select(CRcomparison, 
                            !c(biosense, 
                               stable, 
                               CR_Eco1, # we remove the C:R values as these 
                               CR_Eco2, # are no longer needed: the information
                               CR_MetaEco # is store in column "CR"
                               )) %>% 
  pivot_longer(.,  
               cols = N1rr:PROD_Crr, 
               names_to = "Compartment",
               values_to = "Ratio") %>% 
  dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | 
                                        Compartment == "FLUX_P2rr" | 
                                        Compartment == "FLUX_C1rr" | 
                                        Compartment == "FLUX_C2rr" | 
                                        Compartment == "FLUX_Prr" | 
                                        Compartment == "FLUX_Crr", 
                                      "Nutrient Flux", 
                                      if_else(Compartment == "PROD_P1rr" | 
                                                Compartment == "PROD_C1rr" | 
                                                Compartment == "PROD_P2rr" | 
                                                Compartment == "PROD_C2rr" | 
                                                Compartment == "PROD_Prr" | 
                                                Compartment == "PROD_Crr", 
                                              "Trophic Productivity",
                                                    "Stock"))) %>% 
  dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | 
                                     Compartment == "P1rr" | 
                                     Compartment == "C1rr" | 
                                     Compartment == "FLUX_P1rr" | 
                                     Compartment == "FLUX_C1rr" | 
                                     Compartment == "PROD_P1rr" | 
                                     Compartment == "PROD_C1rr" , 
                                   "Ecosystem 1", 
                                   if_else(Compartment == "N2rr" | 
                                             Compartment == "P2rr" | 
                                             Compartment == "C2rr" | 
                                             Compartment == "FLUX_P2rr" | 
                                             Compartment == "FLUX_C2rr" | 
                                             Compartment == "PROD_P2rr" | 
                                             Compartment == "PROD_C2rr" , 
                                           "Ecosystem 2", 
                                           "Meta-ecosystem"))) %>% 
  dplyr::mutate(., Fertility = factor(Fertility), 
                Function = factor(Function), 
                Scale = factor(Scale), 
                Compartment = factor(Compartment)) 

We now remove all values of the Response Ratio that are \(< 0\) and \(log_{10}\)-transform them before calculating the median.

Show code

CRcomparison_tables <- dplyr::filter(CRcomparison_long, Ratio>=0) %>%
  dplyr::group_by(., CR, Fertility, Function, Scale, Compartment) %>%
  dplyr::mutate(., LRR = log10(Ratio)) %>% 
  dplyr::summarise(.,
                   # min = min(LRR),
                   median = median(LRR), 
                   # max = max(LRR),
                   .groups = "keep")

Finally, we produce summary tables that show the differences in median \(LRR_X\) when considering all stable parameter sets, parameter sets that produce \(C{:}R > 1\), and parameter sets that produce \(C{:}R < 1\). As before, we produce one table per ecosystem function.

Biomass and Nutrient Accumulation

Show code

CRsummStock <- CRcomparison_tables %>% 
  ungroup() %>% 
  filter(., Function == "Stock") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "C1rr", "Consumers" = "C2rr", "Consumers" = "Crr", "Producers" = "P1rr", "Producers" = "P2rr", "Producers" = "Prr", "Nutrients" = "N1rr", "Nutrients" = "N2rr", "Nutrients" = "Nrr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients", "Producers", "Consumers")) %>%
  pivot_wider(., names_from = c(CR, Fertility), 
              values_from = median) %>%
  gt(., groupname_col = "Scale") %>%
  fmt_number(., columns = 4:9, decimals = 2) %>%
  tab_spanner(label = md("**Along-gradient**"), 
              columns = c("all_Along", 
                          "above_1_Along",
                          "below_1_Along")) %>%
  tab_spanner(label = md("**Against-gradient**"), 
                columns = c("all_Against",
                            "above_1_Against",
                            "below_1_Against")) %>%
  cols_label(Compartment = md("**Compartment**"), 
             "all_Against" = md("*All C:R*"), 
             "above_1_Against" = md("*C:R > 1*"), 
             "below_1_Against" = md("*C:R < 1*"), 
             "all_Along" = md("*All C:R*"), 
             "above_1_Along" = md("*C:R > 1*"), 
             "below_1_Along" = md("*C:R < 1*")) %>% 
  cols_hide("Function") %>% 
  tab_footnote(
    footnote = "Includes parameter sets producing C:R > 10.",
    locations = cells_column_labels(
      columns = c("all_Against", "all_Along")
    )
  ) %>% tab_footnote(
    footnote = "Does not include parameter sets producing C:R > 10.",
    locations = cells_column_labels(
      columns = c("above_1_Against", "above_1_Along")
    )
  ) %>% 
  tab_header(title = md("**Median LRR values of biomass and nutrient stocks**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>% 
  tab_style(style = cell_text(align = "left", size = 18), 
            locations = cells_title())

CRsummStock

Median LRR values of biomass and nutrient stocks
LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1.
Compartment	Against-gradient			Along-gradient
	All C:R1	C:R > 12	C:R < 1	All C:R1	C:R > 12	C:R < 1
Ecosystem 1
Consumers	−0.72	−0.71	−0.72	0.26	0.26	0.26
Nutrients	−0.16	−0.55	−0.15	0.12	0.20	0.08
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	0.17	0.15	0.10	−0.27	−0.13	−0.36
Nutrients	0.11	0.12	0.02	−0.14	−0.09	−0.10
Producers	0.04	0.03	0.10	−0.12	−0.12	−0.15
Meta-ecosystem
Consumers	0.08	−0.15	−0.12	−0.09	0.04	−0.09
Nutrients	0.03	−0.18	−0.04	−0.03	0.06	0.00
Producers	0.01	0.02	0.05	−0.02	−0.05	−0.04
1 Includes parameter sets producing C:R > 10.
2 Does not include parameter sets producing C:R > 10.

Show code

gtsave(CRsummStock, filename = "../Results/CR_summary_Stock.rtf")  

Nutrient Flux

Show code

CRsummFlux <- CRcomparison_tables %>% 
  ungroup() %>% 
  filter(., Function == "Nutrient Flux") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "FLUX_C1rr", "Consumers" = "FLUX_C2rr", "Consumers" = "FLUX_Crr", "Producers" = "FLUX_P1rr", "Producers" = "FLUX_P2rr", "Producers" = "FLUX_Prr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
  pivot_wider(., names_from = c(CR, Fertility), 
              values_from = median) %>%
  gt(., groupname_col = "Scale") %>%
  fmt_number(., columns = 4:9, decimals = 2) %>%
  tab_spanner(label = md("**Along-gradient**"), 
              columns = c("all_Along", 
                          "above_1_Along",
                          "below_1_Along")) %>%
  tab_spanner(label = md("**Against-gradient**"), 
                columns = c("all_Against",
                            "above_1_Against",
                            "below_1_Against")) %>%
  cols_label(Compartment = md("**Compartment**"), 
             "all_Against" = md("*All C:R*"), 
             "above_1_Against" = md("*C:R > 1*"), 
             "below_1_Against" = md("*C:R < 1*"), 
             "all_Along" = md("*All C:R*"), 
             "above_1_Along" = md("*C:R > 1*"), 
             "below_1_Along" = md("*C:R < 1*")) %>% 
  cols_hide("Function") %>% 
  tab_footnote(
    footnote = "Includes parameter sets producing C:R > 10.",
    locations = cells_column_labels(
      columns = c("all_Against", "all_Along")
    )
  ) %>% tab_footnote(
    footnote = "Does not include parameter sets producing C:R > 10.",
    locations = cells_column_labels(
      columns = c("above_1_Against", "above_1_Along")
    )
  ) %>% 
  tab_header(title = md("**Median LRR values of nutrient flux**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>% 
  tab_style(style = cell_text(align = "left", size = 18), 
            locations = cells_title())

CRsummFlux

Median LRR values of nutrient flux
LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1.
Compartment	Against-gradient			Along-gradient
	All C:R1	C:R > 12	C:R < 1	All C:R1	C:R > 12	C:R < 1
Ecosystem 1
Consumers	−0.72	−0.71	−0.72	0.26	0.26	0.26
Producers	0.00	0.00	0.00	0.00	0.00	0.00
Ecosystem 2
Consumers	0.17	0.15	0.10	−0.27	−0.13	−0.36
Producers	0.04	0.03	0.10	−0.12	−0.12	−0.15
Meta-ecosystem
Consumers	0.07	−0.01	−0.07	−0.09	0.02	−0.11
Producers	0.01	0.02	0.04	−0.02	−0.05	−0.04
1 Includes parameter sets producing C:R > 10.
2 Does not include parameter sets producing C:R > 10.

Show code

gtsave(CRsummFlux, filename = "../Results/CR_summary_Flux.rtf")

Trophic Productivity

Show code

CRsummProd <- CRcomparison_tables %>% 
  ungroup() %>% 
  filter(., Function == "Trophic Productivity") %>% 
  dplyr::mutate(., Compartment = fct_recode(Compartment, "Consumers" = "PROD_C1rr", "Consumers" = "PROD_C2rr", "Consumers" = "PROD_Crr", "Producers" = "PROD_P1rr", "Producers" = "PROD_P2rr", "Producers" = "PROD_Prr")) %>%
  dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers", "Consumers")) %>%
  pivot_wider(., names_from = c(CR, Fertility), 
              values_from = median) %>%
  gt(., groupname_col = "Scale") %>%
  fmt_number(., columns = 4:9, decimals = 2) %>%
  tab_spanner(label = md("**Along-gradient**"), 
              columns = c("all_Along", 
                          "above_1_Along",
                          "below_1_Along")) %>%
  tab_spanner(label = md("**Against-gradient**"), 
                columns = c("all_Against",
                            "above_1_Against",
                            "below_1_Against")) %>%
  cols_label(Compartment = md("**Compartment**"), 
             "all_Against" = md("*All C:R*"), 
             "above_1_Against" = md("*C:R > 1*"), 
             "below_1_Against" = md("*C:R < 1*"), 
             "all_Along" = md("*All C:R*"), 
             "above_1_Along" = md("*C:R > 1*"), 
             "below_1_Along" = md("*C:R < 1*")) %>% 
  cols_hide("Function") %>% 
  tab_footnote(
    footnote = "Includes parameter sets producing C:R > 10.",
    locations = cells_column_labels(
      columns = c("all_Against", "all_Along")
    )
  ) %>% tab_footnote(
    footnote = "Does not include parameter sets producing C:R > 10.",
    locations = cells_column_labels(
      columns = c("above_1_Against", "above_1_Along")
    )
  ) %>% 
  tab_header(title = md("**Median LRR values of trophic compartment productivity**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>% 
  tab_style(style = cell_text(align = "left", size = 18), 
            locations = cells_title())

CRsummProd

Median LRR values of trophic compartment productivity
LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1.
Compartment	Against-gradient			Along-gradient
	All C:R1	C:R > 12	C:R < 1	All C:R1	C:R > 12	C:R < 1
Ecosystem 1
Consumers	−0.72	−0.71	−0.72	0.26	0.26	0.26
Producers	−0.16	−0.55	−0.15	0.12	0.20	0.08
Ecosystem 2
Consumers	0.22	0.20	0.22	−0.47	−0.28	−0.57
Producers	0.18	0.17	0.16	−0.34	−0.27	−0.34
Meta-ecosystem
Consumers	0.03	−0.07	−0.14	−0.04	0.06	−0.05
Producers	0.03	0.01	0.03	−0.04	−0.02	−0.05
1 Includes parameter sets producing C:R > 10.
2 Does not include parameter sets producing C:R > 10.

Show code

gtsave(CRsummProd, filename = "../Results/CR_summary_Prod.rtf")

This notebook was compiled in 3.7 minutes.